{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# 图嵌在网页里\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ploynomial Curve Fitting "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_toy_data(f, sample_size, std=0.25, noise=True):\n",
    "    x = np.linspace(0, 1, sample_size)\n",
    "    if noise:\n",
    "        y = f(x) + np.random.normal(scale=std, size=x.shape)\n",
    "    else:\n",
    "        y = f(x)\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.sin(2 * np.pi * x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, y_train = create_toy_data(f, 30)\n",
    "x_test, y_test = create_toy_data(f, 100)\n",
    "x_exact, y_exact = create_toy_data(f, 100, noise=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gU1dvG8e9JSIAEQg2d0EvoJKFJqII0la40C1JflSI/ROm9Kh0UEEFFFBSQJkV6b6F3QUIJNYRe0s/7xwCCJCSQ3Z3N7vO5rlxmdmbn3JPgk9kzM+corTVCCCEcn4vZAYQQQtiGFHwhhHASUvCFEMJJSMEXQggnIQVfCCGcRAqzA7xI5syZdd68ec2OIYQQycbevXuva62941pn1wU/b968BAUFmR1DCCGSDaXUufjWSZeOEEI4CSn4QgjhJKTgCyGEk7DrPnwhhOOIiooiJCSE8PBws6M4hFSpUpErVy7c3NwS/R4p+EIImwgJCSFt2rTkzZsXpZTZcZI1rTVhYWGEhISQL1++RL9PunSEEDYRHh5OpkyZpNhbgFKKTJkyvfSnJSn4QgibkWJvOa/ys5QunWTizBnYuROyZIGaNcFF/lQLIV6SRcqGUmqWUuqaUupIPOuVUmqSUuq0UuqQUsrPEu06i1WroEIFWLwYPvsMWrSA2FizUwkhkhtLnSf+ANR9wfp6QKFHXx2Bby3UrlPo0AEWLIDffoO9e42z/aVLzU4lhEhuLFLwtdabgRsv2KQh8JM27ATSK6WyW6JtRxcbC5cvw2uvGcvu7hAQABcvWmDfOpaQOyEcCz3Gvsv72HFhB0evHSXsQRgyE5pwJq89/h/sBR4+fEi1atWIiYnhwoUL1KhRA19fX4oXL87EiRNfqd3IyEiqVq1KdHT0K73/ZdmqDz8ncOGp5ZBHr122UfvJlouL0Z0zZgz06QOnThln9+3avfy+Lt+9zF///MW64HUcvnaYv8P+5kHUgzi3dXNxo5h3MSrkrEDFXBWpXaA2ubxyJfFohLBP27dvT3CbWbNm0aRJE1xdXUmRIgVjx47Fz8+Pu3fv4u/vT+3atSlWrNhLtevu7s7rr7/O/Pnzad269avGTzRbFfy4LifHeQqplOqI0e2Dj4+PNTMlG7/+Co0bw8iRxvLEiVCuXOLee/PhTX4+9DOzD8xm/5X9AGTxzIJ/dn+q56lOkcxFyJg6I6lSpMLd1Z3b4be5cu8Kl+5e4sDVA8w/Op8Z+2YAEOgTyDvF3qFlyZZk9shsjUMVTqL7qu4cuHLAovssk60ME+pOeOE29+/f55133iEkJISYmBj69+/Pu+++S5o0aThy5Aj16tUjMDCQ7du3kzNnTpYsWULq1KkBmDt3Lr/88gsA2bNnJ3t2o5Mibdq0+Pr6cvHiRdzd3alUqRKenp6kT5+e8+fPkyFDBvbv30/Dhg3p06cPtWvXpl+/fty5c4dJkybRqFEjevfu7VAFPwTI/dRyLuBSXBtqrWcAMwACAgKkXwHw8TH67u/cAU9PcHVN+D2nwk4xcutIfj3yK+HR4fhn92fk6yOpU6AOpbOVxkUlrjcvVsdyPPQ4f5z4g/lH59N1VVe+WPsF7cq2o0elHuTLkPiHPoQw26pVq8iRIwd//vknALdv335m/alTp/j111/57rvveOedd1i4cCFt2rQhMjKSM2fOENdw7WfPnmX//v1UqFABLy8vAgMD6dGjB1WqVKF69epMnjwZLy8vBg8ezIABA7h27Rr79+9n6aMLcSVKlGDPnj1WP3bAeGLLEl9AXuBIPOsaACsxzvQrArsTs09/f38tXs65W+d0+yXttetgV+0x3EN3WtZJ772012L7P3jloG67uK12G+KmXQa76LaL2+rLdy9bbP/CcR07dszsCPrkyZM6b968ulevXnrz5s1PXvf09NTBwcG6YMGCT14bNWqUHjp0qNZa64sXL+oiRYo8t7+7d+9qPz8/vXDhwievFSpUSN+4cUNrrXXWrFl1RETEk3VVq1bVfn5++s6dO8/sJ0eOHM+9lhhx/UyBIB1PTbXUbZm/AjuAIkqpEKVUO6VUZ6VU50ebrADOAKeB74CPLdGu+FdkTCTDNw+n8OTC/HToJz4p9wn/dP2HaW9Owy+75e6CLZW1FLMaziK4WzDdKnTj50M/U3hyYcZsG0NEdITF2hHCGgoXLszevXspWbIkvXv3ZsiQIc+sT5ky5ZPvXV1dn1xMTZ069XNPtUZFRdG0aVNat25NkyZNAOPCbnh4OBkyZODChQtkypQJd3d3AA4fPszly5dJmTIladOmfWZfERERpEqVyuLH+1+WukunpdY6u9baTWudS2v9vdZ6mtZ62qP1Wmv9ida6gNa6pNZaZjWxoN0XdxMwI4B+G/rRsGhDTnU5xcR6E8mWJpvV2szplZNxdcZx9OOjVMtbjS/WfoH/DH8OXjlotTaFSKpLly7h4eFBmzZt6NmzJ/v27UvU+zJkyEBMTMyToq+1pl27dvj6+tKjR48n2x07dgxfX18Ajh8//uT7y5cv07p1a5YsWYKnpyerV69+8p6wsDC8vb1fahC0VyXPayZjMbExDN00lErfV+Jm+E2WtljK/Gbz8Ulnu4vdhTIVYlnLZSxvuZywh2GU+64cY7aNISY2xmYZhEisw4cPU758ecqUKcPw4cPp169fot/7xhtvsHXrVgC2bdvGnDlzWL9+PWXKlKFMmTKsWLGCo0ePUqJECcD4VLBv3z6OHTtGkyZNGDt2LL6+vvTv359BgwY92e+GDRuoX7++RY8zXvH19djDl/Thx+/qvau69k+1NYPQrRe21rfDb5sdSYfeD9VN5jfRDELX+qmWvn7/utmRhB2xhz78pNi3b59u06aNxffbuHFjfeLEiVd6ryl9+MK2gi4FUXZ6WTaf28x3b33HnMZz8ErpZXYsMntkZkHzBcx8ayabz22m/MzyHLkW52gbQiQ7ZcuWpUaNGsTEWO7Ta2RkJI0aNaJIkSIW2+eLSMFPZpaeXEq1H6rh7urOrva7aO/X3q5GIFRK0c6vHZs+3MTDqIdUnFmRJSeWmB1LCIv46KOPcE3MfdGJ5O7uzvvvv2+x/SVECn4yMnnXZBrNa0Rx7+LsbLeT0tlKmx0pXhVzVSSoYxDFsxSnyW9N+H7f92ZHEsLpScFPBrTWDNo4iK6rutKwaEM2friRrGmy2qTt+/eNYRx8fKBsWVi7NvHvzZE2B+vfX88bBd6g/bL2jNk2xnpBhRAJkoJv57TW9F3fl8GbBtOubDsWNF+Ah5uHzdrv2NEo+hs3wtCh0LIlHDuW+Pd7unuypMUSWpRowRdrv6D32t4yMJsQJpEJUOyY1prP13zO2B1j6ezfmakNpiZ6SARLWbIELlyADBkgf3545x346y94mTGi3F3d+bnxz3i5ezFq2yjcXN0YUmNIwm8UQliUFHw71nd9X8buGEuX8l2YWHeiKRdn06aF8+eNgg/G9/7+L78fVxdXvn3zW2J0DEM3DyWla0r6Vu1r2bBCiBeSLh07NXb7WEZuHUkn/05WK/Z798LChXD6dPzbDB8Ob74JQ4YYM22dP2+c5b8KF+XC9Den06ZUG/pt6Me4HeNebUdCvIJbt27xzTffvNJ769evz61bt1657TRp0rxwfVKyvQwp+HZo9v7Z9FzTk+bFmjO1/lSrFPvevY0hl3/6CSpVMoZgjstHH8GPP0JEBJQvD1u3QgL/dl/I1cWV2Q1n06xYM/731//49XA8DQthYS8qqgndW79ixQrSp09vjViA7Qq+6U/TvujLGZ+0XX5yuXYZ7KJr/1Rbh0eFW6WNffu0zpVL60cD+unDh7X28tL64UOrNBen8KhwXXV2Ve02xE1vCN5gu4aFaV72SdsNG7QuUEBrd3etq1TR+ty5pLX/7rvv6lSpUunSpUvrnj176g0bNujq1avrli1bal9fX6211g0bNtR+fn66WLFievr06U/emydPHh0aGqqDg4N10aJFdfv27XWxYsV07dq19YMHD55r68yZM7pixYo6ICBA9+vXT3t6emqtjdE1a9asqcuWLatLlCihFy9eHGe2+Lb7r5d90tb0ov6iL2cr+AevHNRpRqTRftP99N2Iu1ZrZ/FirRs0ePa1rFm1vnDBak3G6caDG9p3iq9ONzKdPnL1iG0bFzb3MgX//HmtM2fWeuVKre/f13rYMK3LlNE6NvbV2w8ODtbFixd/srxhwwbt4eGhz5w58+S1sLAwrbXWDx480MWLF9fXrxvDgzxd8F1dXfX+/fu11lo3b95cz5kz57m23nrrLf3jjz9qrbWeMmXKk4IfFRWlb982hkEJDQ3VBQoU0LGxsc9li2+7/5KhFZKpK/eu8OYvb+KV0oulLZaSxj0J/SYJKFUKdu2CA48mHJo/H1KmhGzWG1wzThlSZ2Bl65V4uHnQ4JcGhN4PtW0AYbd27YLKlaFuXfDwMKb3DA6GGy+aOfsVlC9fnnz5/p3EZ9KkSZQuXZqKFSty4cIFTp069dx78uXLR5kyZQDw9/fn7Nmzz22zbds2WrZsCcB777335HWtNX369KFUqVLUqlWLixcvcvXq1efen9jtXpbcpWMHHkY9pOG8hoQ9DGNL2y3k9Mpp1fby5YNp06B6dUiRwphFa/Fi43tby5M+D8taLiNwdiDNf2/OmvfW4OZq/WFihX3LnBn+/hsiI8HdHc6dg6go464xS/L09Hzy/caNG1m7di07duzAw8OD6tWrPzcGPjw/Zv7Dhw/j3Hdc197mzp1LaGgoe/fuxc3Njbx588bZRmK3e1lyhm8yrTWd/+zM7ou7mdtkrkUnK3mRpk0hNBSOHjXOnMqWtUmzcfLP4c/Mt2ay6dwmPlv9mXlBhN2oWhVKlIDAQOjSxVgePdoo/q8qbdq03L17N971t2/fJkOGDHh4eHDixAl27tz5ym1VrlyZefPmAUbxfrqNLFmy4ObmxoYNGzh37lyc2eLbLqmk4Jts+t7p/HTwJwZWG0ijoo1s2rabG2TNCi528K+gdanW9KzUk6l7pvLd3u/MjiNM5uIC8+bBF19AgQLGXWSffpq0fWbKlInKlStTokQJPv/88+fW161bl+joaEqVKkX//v2pWLHiK7c1ceJEpk6dSrly5Z6ZN7d169YEBQUREBDA3LlzKVq0aJzZ4tsuqZS248fcAwICdFCQ406OtStkF1VmV6FW/losb7Xc5k/R2puY2Bjq/1KfjWc3sv2j7fjneIUnvITdenoGKGEZcf1MlVJ7tdYBcW3v3BXGSsLDYc0aWL3aGIcmLqH3Q2n2ezNyeuXk5yY/O32xB+Me/blN5pLFMwvvLHiH2+G3E36TECLR5KKthd24ATVqGHcWpEhh9JNv3PjsHTCxOpb3F79P6P1QdrTbQcbUGU3La28ye2RmfrP5VPuhGu2WtuP35r/b1Xj/QiRnclppYUOHwmuvwfbtsGULvP02/HfazAk7J7Dq9CrGvjGWstlNvFpqp17L/RojXx/JwuMLmbJ7itlxhAXZcxdycvMqP0sp+C8hLAw2bYITJ+Lf5tw5qFkTHp+U1qxpvPbY3kt7+XLtlzQq2oiPy31s3cDJ2P8q/Y+3Cr9FzzU9OXjloNlxhAWkSpWKsLAwKfoWoLUmLCyMVKlSvdT7pEsnkbZsgWbNoFAh+OcfeP994zax/woIgFmzjAHHXF1h5sx/R5e8G3GXFgtbkDVNVr5/+3un7ar4+2/jjot//jFuB5061bhb6GlKKWY1nEXJb0vSelFr9nTYQ2q31OYEFhaRK1cuQkJCCA2VB+wsIVWqVOTKlevl3hTfI7j28GVPQyv4+Gi9YoXx/c2bxhgfGzc+v11kpNYtW2qdNq3W6dJp3bCh1o+H2vho8UfaZbCL3nR2k+2C25lbt7TOnVvriRO1PnlS688/1zogQOuYmLi3X3VqlWYQuuuKrrYNKkQyhQytkDSRkXDxovGYN0D69FClCsTx1DVubvDLL8bDTKdPG0+wpk4Ny04uY9aBWXxR+Quq5qlqtawPH8KKFbB0Kdy2w5tcgoIgb17o2hUKFzY+JV28aEyyEpc6BevQtXxXJu2exOrTq22aVQhHIwU/EdzdjYc/fvnFWL54Edatg5Il439PpkzG4+EA1x9cp8OyDpTKWoqB1QZaLefNm8ZQxyNGwJQpUKaMMX69PfH0hGvXIDraWL5zx7h19akn3J8zqtYoinsX58MlH3LjoYUHUxHCiUjBT6T5840x5AsXNqb369YNKlRI+H1aazov78yNhzeY03gOKVP8Ow7H0aPwxhvg6wtt20IS5lcAYNQoKFfOuN7w11/GPr/8Mmn7tLTy5aFgQWjQAL76CmrVMq6HPP7jGJfUbqmZ03gO1x9cp9uqbrYLK4SDkYKfSGXKGBcblywxLjb+73+Je9+vR35l4fGFDK0xlFJZSz15PTQUatc2JiH5/Xfjnv1XnUnqsfPnjTFHHl8Lrlo1/q4Sa5g1C4oWNQZn698f4ppTwsUF/vgD3noLLl82unYmTUp432Wzl6VPYB9+PvQzS08utXx4IZyADK1gRdfuX6PY1GIUzFiQbR9tw9XF9cm6RYuMArl8ubEcE2NcGwgJgXTpXq29CROMYrp8uTHccevW4OMDY8da4GASsHSp8ann11+N/O3aGUW9d2/LtREZE0n578pz9f5Vjn58VB5YEyIOMrSCSbqu7MrdyLvMajjrmWIPxpO4oaHw+O/tzZtGv/ZTI6++tC5doEgR4xbHjBmN/Q0dmoQDeAlLl0KvXlCxotFFNXq08WnIktxd3fmh0Q9cf3Cdriu7WnbnQjgBKfhWsuTEEuYfnU//qv0p5l3sufU1axrdOM2bw/jxRl92t27wks9RPMPVFWbMMC6KXrxonO17eCThIF6Cl9ezF4jPnzdes7Qy2crQJ7APcw/PZeWplU9ej42FgweNiTPiGZ5cCKcnXTpWcCv8FsWmFsPb05ugDkHxTujx4IFxN01IiHEBuFWrf/vfk5tz54w7hBo1Mrp0Zs40uq2qVLF8WxHREZSdXpYHUQ848vER3ElD06bGRfB06f4dvO5ln0kRwhG8qEtHnrS1gi/WfMHV+1dZ1nLZC2dv8vAwukEcQZ48sHs3/Pij8dzCunXGVIrWkDJFSma8NYMqs6swYMMAch8bR3Q0nDxpPAcxaBB07w4LFlinfSGSKyn4Frb9wnZm7JtBj4o9nG4891y5oG9f27QV6BNIZ//OTNw1kUahLWnQoBxuj/62NmoECxfaJocQyYn04VtQVEwUnZZ3IrdXbgbXGGx2HIc3qtYosnpmZXfWjixYFE14uHERfO5cY3o8IcSzpOBb0Pid4zly7QiT600mjXsas+M4vHSp0jG53mRCYg5wt+hU8uQxHupat864RVUI8Szp0rGQs7fOMmjjIBoWaUjDog3NjuM0mvg2oW7Bumxz78/aDe+QVmWnUCHjDighxLPkDN9Cuq3qhotyYVK9RDw2KixGKcXkepOJjIlkwvH/4esrxV6I+Fik4Cul6iqlTiqlTiulnhu9RSlVXSl1Wyl14NHXAEu0ay9WnFrB0pNL6V+1Pz7pfMyO43QKZizIl4Ff8uuRX1l3Zp3ZcYSwW0m+D18p5Qr8DdQGQoA9QEut9bGntqkO9NRav/ky+04O9+FHREdQ4tsSuCpXDv3fIdxd3c2O5JQeRj2k5LclSeGSQn4PwqlZe2iF8sBprfUZrXUkMA9wmk7sr7d/zekbp5lcb7IUGROldkvNpHqTOBl2ksm7Jse73YED4OcHadMacw/HNaeBEI7KEgU/J/D0mIwhj177r0pKqYNKqZVKqeLx7Uwp1VEpFaSUCrL3qdDO3z7P8C3DaerblNoFapsdx+nVL1SfBoUaMHjTYK7cu/Lc+tu3jWGZP/vMGEW0ZUuoV894UEwIZ2CJgh/XYAD/7SfaB+TRWpcGJgOL49uZ1nqG1jpAax3g7e1tgXjW0/OvngCMqzPO5CTisfF1xhMeHU7vdc8P03nwoPFE8HvvGSOTduliDGXxzz8mBBXCBJYo+CFA7qeWcwGXnt5Aa31Ha33v0fcrADel1AumvLB/m89t5vdjv/Nl4JdyodaOFMpUiB6VevDDgR/YFbLrmXWZMhln9vfvG8thYXD9OmTIYEJQIUxgiYK/ByiklMqnlHIHWgDPzFChlMqmlDEsmFKq/KN2wyzQtiliYmPovqo7ub1y0/O1nmbHEf/Rt0pfsqfJTpeVXYjVsU9eL1bM6MIJDISePaFyZfj4Y8iWzcSwQthQkgu+1joa+BRYDRwHftNaH1VKdVZKdX60WTPgiFLqIDAJaKHteZjOBMw+MJv9V/YzpvYYPNxsNP6wSLS0KdMyutZo9lzaw8+Hfn7yulIwfboxG5e3tzEs9fDhJgYVwsZkeOSXdCfiDoUmF6JQxkJsabsFlVzHM3ZwsTqWCjMrcOnuJf7+9G883V8wS7oQDkRmvLKg4ZuHE3o/lIl1J0qxt2MuyoXxdcZz6e4lvtr+ldlxhLALUvBfQvDNYCbsmsD7pd93uqGPk6NAn0CaF2vOmG1jCLkTYnYcIUwnBf8l9F7XG1flyvCa0vGbXIyuNZoYHUOfdX3MjiKE6aTgJ9LOkJ3MPzqfnq/1JKdXXM+VCXuUL0M+Pqv4GXMOzSHokn1dDxLC1qTgJ4LWmh6re5AtTTZ6VXaQOQmdSO/A3mT2yMznaz7Hnm9SEMLapOAnwoJjC9gRsoNhNYbJxCbJULpU6RhYbSAbz27kz1N/mh1HCNPIbZkJiIyJxHeqL55unuzvtB9XF1dT84hXExUTRfFvij8ZTTOFiwyaLxyT3JaZBNODpnPm5hlG1xotxT4Zc3N1Y1StURy/fpzZ+2ebHUcIU0jBf4E7EXcYsnkINfPVpG7BumbHEUnUuGhjKueuzICNA7gXec/sOELYnBT8FxizbQzXH1xnTK0x8pCVA1BK8VXtr7hy7woTdsos58L5OH3BDw42xkWvVg369IGICOP1i3cuMm7HOFqWaCkPWTmQSrkr0bhoY8ZsG0Poffueb0EIS3Pqgh8WZhT6kiVh4EA4fBjatTPWDdo4iOjYaHnIygGNeH0E96PuM3yL/G6Fc3Hqgr92LZQpY5zZ16wJv/0Gv/8Ohy6dYNaBWfxfwP+RL0M+s2M6tdBQWLECduwAS91QVjRzUdqVbcc3e74h+GawZXYqRDLg1AXf1RXCw/9dftydM3hLfzzcPOhbta85wQQAe/YYn77Gj4cPPoB334WYGMvse2C1gaRwSUH/Df0ts0MhkgGnLvh16hgzIHXpAnPnGvOdNvp4D4tOLKBnpZ5k8cxidkSn1rEjTJgAa9YY3W0XLsD8+ZbZd06vnHSv2J25h+dy4MoBy+xUCDvn1AU/bVrYsgXc3GDZMmjRAsLK9sbbw5selXqYHc/pnT0LtWoZ36dMaVxvCbZgD0yvyr3IkCoDfdfLJzlhvthYy3VbxsepCz5A5swwbhzMmwdF669hXfA6+lXtR9qUac2O5vT8/GDaNON/gqtXYdEi8LfgDVPpU6Xni8pfsOLUCrae32q5HQvxEqKioHNn8PSENGmMa4rWKvxOX/Af01rTZ30f8qTLQyf/TmbHEcDs2caF9Bw5oEABeO89qGvh59+6VOhC9jTZ+XLtlzKwmjDFiBFw+jRcumT896+/YMYM67QlBf+RP078QdClIAZXH0zKFCnNjiMAHx84cMC4eHvpkjEXraV5uHkwoNoAtl3YxopTKyzfgBAJ2LABeveGDBkge3bo3t14zRqk4AMxsTH0W98P38y+tCnVxuw44ikuLpArF3h5Wa+NdmXbUSBDAfqs70OsjrVeQ0LEIUsW48TmsYMHjdesQYYMBH4+9DPHrx9n4TsLZYA0J+Tm6saQGkNovag1vx39jRYlWpgdSTiRoUOhenXYu9e4TXz/fti+3TptOf3wyBHRERSZUoTMHpnZ02GPjJnjpGJ1LGWmlSE8OpxjnxyT4ZOFTV2+DH/+aTwb1LAhZMz46vuS4ZFf4Lt933Hu9jlGvD5Cir0Tc1EuDK0xlFM3TvHjgR/NjiOcTPbs0L49tG2btGKfEKcu+A+iHjB8y3Cq5alG7fy1zY4jTPZ2kbcpn7M8gzcNJiI6ItHv27wZ6tUznhOYNMn691IL8aqcuuBP2T2FK/euMLzmcDm7FyilGF5zOBfuXGD63umJes/+/dCsGbRqBf36waxZ8PXXVg4qxCty2oJ/O/w2o7eNpl7BelT2qWx2HGEnXs/3OtXzVmf4luHcj7yf4Pa//gqffGI8I1C7NsycCT9Kj5CwU05b8CfsnMCNhzcYWmOo2VGEHXl8ln/t/jWm7J6S4PZubnD/qb8L9+9DCrneKx558ACOHYObN81OYnDKgh/2IIyxO8bSxLeJTG4invNa7teoX6g+o7eN5nb47Rdu+9FHxhn9kCHw3Xfw/vvQQ4ZhEsC2bZA/PzRuDPnyGZ/+zOaUBf+r7V9xL/IeQ6oPMTuKsFNDqg/hZvjNBKdCLFDAuGh75Qps3QqTJxtFXzi36Gjj2s6sWXDyJAQFQd++xvdmcroPn1fuXWHSrkm0KtmK4lmKmx1H2Cn/HP40LtqEr7ePo0WBTymSO1O82xYpAt98Y8Nwwu6FhhpzN9SvbywXLAgVKhjdO0WKmJfL6c7wR20dRWRMJAOrDTQ7irBjFy7A4UmDuRd5l1Kdv+Lzz+V2S5F4mTMbBX/HDmP5yhXjLL9gQXNzOWTB1xpGjzb6zfLnh6++Ml4LuRPCtKBpfFD6AwplKmR2TGHHOnSA9+qUoFXJlri+NonlG6+yaJHZqURy4eYGc+bA229DlSpQqhR07WrM4GYmh+zS+e47YwarpUuN5VatjJHo9mYbTqyOpX81mdZOvNiBA/D993A/1UDmHZ1HugajOHBgPE2bmp1MJBf168ORI3D8uDEAoNln9+CgZ/jLlsHgwcZf05IlYeBAmLc6mO/3f097v/bkTZ/X7IjCzhUsaEyeXjhTYdqU+IDdsd+SIU+I2bFEMpM1qzEwmj0Ue3DQgu/lZUyP99i5c3DWZyguyoW+VWQ6O5GwadOMk4aqVWHD4P7gEsOpbCPMjiVEkjhkl07v3lCzpuOjjR0AABy8SURBVFHotYY5K05x572f6BLQhZxeOc2OJ5KBEiWMj+P79kG6dPn47nJ7vt8/ky8Ce8knRJFsOeQZfokSxtXxzJmNiQSq9jdmsfoy8Euzo4lkJH1648TB3x/6Ve1rjKi5SZ7MFsmXQxZ8MB6I6dcPGnc8xtLgX/i03KdkTZPV7FgimcrllYtO/p348eCPnAo7ZXYcIV6JRQq+UqquUuqkUuq0Uuq502hlmPRo/SGllJ8l2k2MQRsH4enuyeeVP7dVk8JB9a7SG3dXd4Zslie0RfKU5IKvlHIFpgL1gGJAS6VUsf9sVg8o9OirI/BtUttNjANXDvD7sd/pXqE7mT0y26JJ4cCypcnGJ+U+Ye6huRwPPW52HOFg7kbcxdozEFriDL88cFprfUZrHQnMAxr+Z5uGwE/asBNIr5TKboG2X2jgxoGkS5mOHpVkNCthGb0q98LT3ZNBmwaZHUU4mMbzG/Pugnet2oYlCn5O4MJTyyGPXnvZbQBQSnVUSgUppYJCQ0NfOZTWmjcLvcmI10eQIXWGV96PEE/z9vSmW4Vu/Hb0Nw5dPWR2HOEgNp3dxLrgdVTKVcmq7Vii4Mc1VdR/P5ckZhvjRa1naK0DtNYB3t7erx5KKTr4d+Djch+/8j6EiMv/Kv2PdCnTMXCjjMckkk5rTf8N/cmRNgedAzpbtS1LFPwQIPdTy7mAS6+wjRDJQobUGehRqQeLTyxm76W9ZscRydzaM2vZcn4LfQL7kNottVXbskTB3wMUUkrlU0q5Ay2Apf/ZZinw/qO7dSoCt7XWly3QthCm6F6xOxlTZ2TAxgFmRxHJ2OOz+9xeuWnv197q7SW54Guto4FPgdXAceA3rfVRpVRnpdTjzycrgDPAaeA7QPpZRLLmldKLz1/7nBWnVrDjwg6z44hkasWpFey6uIv+VfuTMkVKq7enrH0bUFIEBATooKAgs2MIEad7kffIPzE/pbOVZs17a8yOI5IZrTUB3wVwK/wWJz45gZurm0X2q5Taq7UOiGudwz5pK4S1pXFPQ+/A3qw9s5ZNZzeZHUckM3+c+IN9l/cxoOoAixX7hEjBFyIJOgd0JkfaHPTf0N/qD80IxxETG8OADQMomrkobUq1sVm7UvCFSILUbqnpE9iHLee3sOaMdOuIxPnt6G8cDT3KoGqDcHVxtVm7UvCFSKL2fu3xSecjZ/kiUaJjoxm4cSAls5SkefHmNm1bCr4QSZQyRUoGVB3A7ou7Wf73crPjCDv386GfOXXjFENqDMFF2bYES8EXwgLeL/0+BTMWpP+G/sTqWLPjCDsVGRPJ4E2D8c/uT8Mi/x1yzPqk4AthAW6ubgyqNoiDVw+y8NhCs+MIO/X9vu85e+ssw2oOQ6m4RpyxLin4QlhIixItKOZdjAEbBxATG2N2HGFnHkY9ZNiWYQT6BFKnQB1TMkjBF8JCXF1cGVJ9CCeun2Du4blmxxF25tugb7l09xLDaw435ewepOALYVFNfJvgl92PQRsHERkTaXYcYSfuRtxl5NaR1M5fm6p5qpqWQwq+EBaklGJYjWEE3wrm+33fmx1H2IlJuyZx/cF1htUcZmoOKfhCWFjdgnUJ9Alk6OahPIx6aHYcYbIbD2/w1faveLvI25TPWd7ULFLwhbAwpRTDaw7n8r3LTN0z1ew4wmRfbfuKOxF3GFbD3LN7kIIvhFVUzVOVNwq8waito7gTccfsOMKCYmJg2TKYPRv+/vvF2165d4WJuybSqmQrSmYtaZuALyAFXwgrGV5zOGEPwxi3Y5zZUYSFREdDw4YweDCsXw+VK8PKlfFvP2zzMKJioxhUfZDNMr6IFHwhrCQgRwBNfZsydsdYQu+Hmh1HWMCCBXDzJuzcCXPmwO+/w8fxTOcUfDOYGXtn0K5sOwpmLGjboPGQgi+EFQ2tMZQHUQ8YuXUkWsOKFTBuHKxaZXYy8SquXAE/P0iRwlguVw4uxTM79+BNg3F1caV/1f62C5gAKfhCWJGvty8flv6Qb/Z8Q4ee5/n8czh3Drp3hy+/NDudeFmVKsHChXD8OMTGwvDhUKXK89sdvXaUOYfm8Em5T8jpldP2QeMhBV8IKxtYfSCxWjM3ZDDbt8PEibBjB8ycCefPm51OvIwKFYwiX6ECpEoFW7YYXTv/1Wd9nyczotkTKfhCWJlPOh+a5v6EcN8fuBx1AoAMGSB7dggLMzmceGlt28Lt23DnjlHws2d/dv32C9tZenIpvV7rRSaPTOaEjIcUfCFsYES93qgoT96b3Zf792HWLKNoFClidjLxKpQyzvD/S2vNl2u/JFuabHSv2N32wRIgBV8IG8iX1Zv/K/U5QQ8Wkb7ETqZMMS7geniYnUxY0srTK9lyfgsDqg7A093T7DjPUfY8JVtAQIAOCgoyO4YQFnEv8h4FJxWkcKbCbPpwk2kjJgrriImNwW+GH/cj73P8k+O4ubqZkkMptVdrHRDXOjnDF8JG0rinYWC1gWw5v4UVp1aYHUdY2NzDczl09RDDag4zrdgnRM7whbChqJgoin1TjFQpUnGg0wFcXVzNjiQsIDw6nMKTC5M1TVZ2td9l87lqnyZn+ELYCTdXN0bUHMGRa0eYcyiO+/lEsjR512Qu3LnAmFpjTC32CbHfZEI4qGbFmlE+Z3n6b+jPg6gHZscRSXTj4Q1GbB1B/UL1qZGvhtlxXkgKvhA2ppTi69pfE3InhAk7J5gdRyTRiC0juB1+m1GvjzI7SoKk4Athgip5qtCwSENGbR3FtfvXzI4jXlHwzWAm757MB2U+sIvhjxMiBV8Ik4yuNZoHUQ8YsmmI2VHEK+q9rjeuytUuJjdJDCn4QpikSOYidPTvyPS90zl5/aTZccRL2hmyk/lH59PztZ52NUDai0jBF8JEg6oPInWK1PRa28vsKOIlaK3psboH2dJko1fl5PO7k4IvhImyeGahT5U+LD25lPXB682OIxJpwbEF7AjZwdAaQ0njnsbsOIkmD14JYbLw6HB8p/rildKLfR33ycNYdi4iOgLfqb54unva5cNz8uCVEHYsVYpUjK41mkNXDzH7wGyz44gETNg5geBbwYx9Y6zdFfuESMEXwg40L9acyrkr0299P+5G3DU7jojH5buXGbZlGG8Vfos3CrxhdpyXJgVfCDuglGJ8nfFcvX+V4VuGmx3HKc2dC3Xrwltvwbp1cW/Td31fIqIjGPvGWNuGs5AkFXylVEal1Bql1KlH/80Qz3ZnlVKHlVIHlFLSKS9EHMrlLMeHZT5k/M7xnAo7ZXYcp/LTTzBgAHTuDO+8Ay1bwrZtz24TdCmIHw78QPeK3SmUqZA5QZMoqWf4XwLrtNaFgHWPluNTQ2tdJr6LCUIIGPn6SFK6pqTHXz3MjuJUZs6EqVOhUSN47z3o0wd+/PHf9Vpruq/qjrenN/2q9jMvaBIlteA3BB7/WH4EGiVxf0I4tWxpsjGw2kCW/71cxsy3IVdXiIz8dzkyElyeqo5zD89l24VtjKg5Aq+UXrYPaCFJui1TKXVLa53+qeWbWuvnunWUUsHATUAD07XWM16wz45ARwAfHx//c+fOvXI+IZKjyJhISn1bilgdy5GPj+Du6m52JLuzbRv07Qs3bhj97sOGgXsSfkwLF0K3bjB0qDE5+bBhsHo1+PnBnYg7FJlSBJ90Puxot8Ouhz+GJN6WqZRaq5Q6EsdXw5fIUFlr7QfUAz5RSlWNb0Ot9QytdYDWOsDb2/slmhDCMbi7ujOh7gRO3TjF+B3jzY5jd06cMLpeOnWCH36AQ4egexLnC2/aFKZPh7/+gn37jPmG/fyMdYM3DubqvatMqTfF7ot9QpJ6hn8SqK61vqyUyg5s1FoXSeA9g4B7WuuvE9q/PHglnFmjeY1Yc2YNxz85jk86H7Pj2I2vv4Zz52DyZGP50iUoWRLCwizf1tFrRyk9rTQflf2IGW/F2zFhV6z54NVS4INH338ALImjcU+lVNrH3wNvAEeS2K4QDm9i3Ylorfls9WdmR7ErqVI9W9yvX4fUqS3fjtaaLiu74JXSixGvj7B8AyZIasEfBdRWSp0Caj9aRimVQyn1+IpTVmCrUuogsBv4U2u9KontCuHw8qTPw4BqA1h0fJFcwH1Kq1awaxd8/DFMnGh07/Tta/l2fjn8CxvObmB4zeFk9shs+QZMIGPpCGHHImMiKT2tNJExkRz5vyOkdrPCqWwydO2a0aVz4wbUqQNvv23Z/d98eJOiU4uSN31etn+0PVkNoSBj6QiRTLm7ujO1/lTO3DzDiC2O0a1gCVmyGHfUTJ1q+WIP8OXaLwl7EMb0N6cnq2KfECn4Qti5mvlq0qZUG0ZvG83Ra0fNjuPwtl/Yzox9M+hWoRtlspUxO45FScEXIhkY98Y4vFJ60WFZB2J1rNlxHFZUTBSdl3cmt1duBtcYbHYci5OCL0Qy4O3pzfg649kRsoNv93xrdhyHNXrbaA5fO8yU+lOS1cQmiSUFXwg7FRsLT99T0aZUG2rnr03vdb0JuRNiXjAHdSz0GEM3D+Xd4u/ydhErXBiwA1LwhbAz9+5B8+bGveUZM8KUKcbrSimmvTmN6NhoOi/vjD3fYZfcxMTG0G5pO9K6p2VSvUlmx7EaKfhC2JmuXSFlSrh5E3bvhrFjYdWjJ1fyZ8jPiNdH8OepP/np4E/mBk0m1q83foaLFj37ielpk3dPZmfITibWnUgWzyy2DWhDUvCFsDMbNsCgQeDhAYUKQYcOxmuPdSnfhcq5K9NtVTcu3rloWs7kYORIaN8eLlwwBkRr2/b5on/6xmn6ru9Lg0INaFWylTlBbUQKvhB2JksWOHDA+F5rOHjQeO0xVxdXZjecTWRMJJ2Wd5KunXjcvg0jRhgja06YYPx382Z4+lnOmNgYPlj8AW4ubkx7cxpKKfMC24AUfCHszNix8H//Bx99ZDxFeuqUMTLk0wplKvSka+fHgz/GvSMnd+sWeHlB9uzGcurUkD//s+PwfL39a7Zf2M7U+lPJ5ZXLnKA2JAVfCDsTGAg7d0KFCvDBB7B1K6SJ4w7BrhW6UjVPVbqu7ErwzWDbB7VzuXIZBX/sWONC+KJFcPgw+Psb6w9dPUT/Df1pVqyZw3flPCZj6QiRjJ27dY7S00pTPEtxNn24iRQuKcyOZFfOnDGmLNy3zzi7//57qFgRIqIjKD+zPFfvXeXIx0ccZnA0kLF0hHBYedLn4ZsG37D9wnZGbhlpdhy7kz+/0Xf/8CEcPWoUezDGyjl09RAz357pUMU+IVLwhUjmWpVsRauSrRi8aTC7QnaZHcfu/fn3n0zYNYEu5bvwZuE3zY5jU1LwhXAAjy86tljYgpsPb5odx25dunuJD5d8SOmspRlTe4zZcWxOCr4QDiB9qvTMazaPkDshtF3SNt5bNaOjbRzMjsTExtBmURseRD1gXrN5pEqRyuxINicFXwgHUTFXRcbUGsOSk0uYsHPCM+v274eiRY0neAsWNGaMcjaDNw1mw9kNTK43maKZi5odxxRS8IVwIN0rdqdR0Ub0WtuLHRd2APDgAbz1FgwcCFFR8NVXxrSAt2+bHNaGlv+9nKGbh9K2TFvalmlrdhzTSMEXwoEopZjdcDa5vXLT7PdmXL57mVOnIH16aNkSXFygcWPInRuOHTM7rW2cvnGaNova4Jfdj6n1pzr807QvIgVfCAeTPlV6FrdYzK3wWzT9rSleGSO4dMmYBxaMeWDPnoWsWU2NaRMPoh7Q9LemuLq4svCdhU4/J7AUfCEcUKmspfip0U/sCNnB8P2f8NlnmgoVjIHEypc3/ps/v9kprStWx/LeH+9x+Oph5jaZS970ec2OZDop+EI4qKbFmtKvSj++3/896d6YzJw5UK4czJxpDCrm6Pqu68ui44sYV2ccdQvWNTuOXZDnsIVwYINrDOZI6BG6r+rOond96NSpkdmRbOKHAz8watsoOvl3oluFbmbHsRtyhi+EA3NRLsxtMpfyOcvTcmHLJ3fuOLJ1Z9bRcVlHauWvxeR6k536Iu1/ScEXwsF5uHmwrOUycnnl4q1f3+LvsL/NjmQ1ey7uodH8RhTJXITfmv2Gm6ub2ZHsihR8IZyAt6c3K1uvRCnFG3Pe4Nytc2ZHsrjjocepN7ce3h7e/NXmLzKkzmB2JLsjBV8IJ1EwY0FWt1nNrfBbvP7T6w41PWLwzWDe+PkNUrikYM17a8ieNrvZkeySFHwhnIhfdj9Wt1nNtfvXeP2n17ly74rZkZLs9I3TVPuhGvcj77O6zWoKZCxgdiS7JQVfCCdTIVcFVrReQcidEKr/UJ0Lty+YHemVnbx+kmo/VONh9EM2fLCB0tlKmx3JrknBF8IJBfoEsrL1Si7fu0zlWZU5cf3Ek3VBQdC/P4wcCVevvng/V67A338bY/TY2qGrh6j2QzWiY6Ol2CeSFHwhnFSVPFXY9OEmImMiCZwVyJ6Le1i9GurXN9YHBxsPal2+/Px7tYb//c8YgbNuXShdGs7Z8Drw6tOrCZwVSAqXFGz8YCMlspSwXePJmBR8IZxYmWxl2PrRVtKmTEu1H6rx6bR5zJwJQ4fCjBnw9tswbdrz7/vjD1izxvij8Hje2HbtbJN55r6ZNPilAfkz5GdX+134evvapmEHIAVfCCdXMGNBdrbbiX8Of06Xacnie72JiY0BwMcH7t59/j2HDkHDhpDh0Z2PH35ovGZNEdERdF3ZlQ7LOlC7QG22tN1CTq+c1m3UwUjBF0KQNU1W1r2/jjIxHZl9ahTVvqvPojWXmDgR3oxj2teCBWHtWggPN5aXLzdes5bgm8EEzg5k8u7JfFbxM5a2WEralGmt16CDkrF0hBAAuLu6s7v/dN4c6M+akO7sPFeSLoNmULNm0+e2bdUKVq40+vBz5IALF2DVKstn0loz9/BcuqzsgtaaP979g0ZFnWM8IGtQ8c19aQ8CAgJ0UFCQ2TGEcDonrp/gvT/eI+hSEK1KtuKr2l+RI22OZ7bRGg4eNGbOKlMG0qWLf38nTxoTrhQoAKVKJS5D8M1g/u/P/2P1P6upmKsic5vMJX8GBx/T2QKUUnu11gFxrZMuHSHEc4pmLsr2j7YzsNpAFhxbQOHJhRm5ZSTh0eFPtlHKKPTVqr242H//PVStCrNnQ716xu2eL3In4g6DNg6ixLcl2HZhG5PrTWZr261S7C1AzvCFEC/0z41/6LmmJ4tPLCa3V256VOpBe7/2pHFPk+B7b92CvHmNe/sLFjTu6y9VCrZte77P/27EXaYFTWPUtlHceHiDpr5NGV9nPLnT5bbOgTkoq53hK6WaK6WOKqVilVJxNvBou7pKqZNKqdNKqS+T0qYQwrYKZCzAH+/+wdr31pIvQz4+W/0ZPuN96LWmF/sv7+dFJ42XL0OWLP8W96xZoUgRo88fjD76fZf30WlZJ3KMy0Gvtb0on7M8QR2CWPDOAin2FpakM3yllC8QC0wHemqtnzsdV0q5An8DtYEQYA/QUmud4BTKcoYvhP3ZGbKTr7Z/xdKTS4mOjaZQxkI0KtqISrkqUSFXhWf6+h8+NKZSnDkTGjSAHTs0DVpeYty8IHaHrWb1P6s5c/MMqVOk5t0S79LZvzMVclUw8eiSvxed4VukS0cptZH4C34lYJDWus6j5d4AWusEevKk4Athz8IehLHo+CLmH53P5nObiYo1xlfw9vAme9rsZEuTDa+UXly+FsGe/RHEpgwjOv1JcL8HgKebJzXz1aR+ofq0KNGC9KnSm3k4DuNFBd8Wt2XmBJ4enSkEiPdPuFKqI9ARwMfHx7rJhBCvLJNHJjr4d6CDfwfCo8PZf3k/uy7u4ui1o1y9f5Ur965w/vZ5UqVIRZmKKfFQmSmRrTK+WYpQIksJKuaqiLuru9mH4VQSLPhKqbVAtjhW9dVaL0lEG3HNLxbvxwqt9QxgBhhn+InYvxBO6+JFmDrVeBr27behdm1zcqRKkYpKuStRKXclcwKIREnwoq3WupbWukQcX4kp9mCc0T995SUXcOlVwgoh/nX5MlSsaPST58sHbdvCL7+YnUrYM1vch78HKKSUyqeUcgdaAEtt0K4QDm3WLGPYg/HjoUcPo9gPH252qlcXHQ1ffGGM31O4sHHfvrCsJPXhK6UaA5MBb+BPpdQBrXUdpVQOYKbWur7WOlop9SmwGnAFZmmtjyY5uRBO7uFDyJz532Vvb+O15GrYMNi5E9atg5s3oXlz45bOBg3MTuY4klTwtdZ/AH/E8foloP5TyyuAFUlpSwjxrMaNjSdX/fwgd27jLP/dd81O9eqWLoVvv4VChYzlnj2NQdmk4FuODJ4mRDLl7w9z58LgwcZF20aNjJmqkqt06eDsWajw6B6+4GDw8jI1ksORgi9EMla7tnl35ljaoEHQrBns2wc3bhijb+7caXYqxyKDpwkh7EK1akb/vaen0a2zZw/klPlNLErO8IUQdqNUqcQPnyxenpzhCyGEk5CCL4QQTkIKvhBCOAkp+EII4SSk4AshhJOQgi+EEE5CCr4QQjgJKfhCCOEkLDLFobUopUKBc0ncTWbgugXiJBdyvI5NjtexWeJ482itveNaYdcF3xKUUkHxze/oiOR4HZscr2Oz9vFKl44QQjgJKfhCCOEknKHgzzA7gI3J8To2OV7HZtXjdfg+fCGEEAZnOMMXQgiBFHwhhHAaDlPwlVJ1lVInlVKnlVJfxrFeKaUmPVp/SCnlZ0ZOS0nE8bZ+dJyHlFLblVKlzchpKQkd71PblVNKxSilmtkyn6Ul5niVUtWVUgeUUkeVUptsndGSEvHvOZ1SaplS6uCj421rRk5LUErNUkpdU0odiWe99WqV1jrZfwGuwD9AfsAdOAgU+8829YGVgAIqArvMzm3l430NyPDo+3qOfrxPbbceWAE0Mzu3lX+/6YFjgM+j5Sxm57by8fYBRj/63hu4Abibnf0Vj7cq4AcciWe91WqVo5zhlwdOa63PaK0jgXlAw/9s0xD4SRt2AumVUtltHdRCEjxerfV2rfXNR4s7gVw2zmhJifn9AnQBFgLXbBnOChJzvK2ARVrr8wBa6+R8zIk5Xg2kVUopIA1GwY+2bUzL0FpvxsgfH6vVKkcp+DmBC08thzx67WW3SS5e9ljaYZwxJFcJHq9SKifQGJhmw1zWkpjfb2Egg1Jqo1Jqr1LqfZuls7zEHO8UwBe4BBwGummtY20Tz+asVqscZRJzFcdr/73fNDHbJBeJPhalVA2Mgh9o1UTWlZjjnQB8obWOMU4Ck7XEHG8KwB94HUgN7FBK7dRa/23tcFaQmOOtAxwAagIFgDVKqS1a6zvWDmcCq9UqRyn4IUDup5ZzYZwJvOw2yUWijkUpVQqYCdTTWofZKJs1JOZ4A4B5j4p9ZqC+Uipaa73YNhEtKrH/nq9rre8D95VSm4HSQHIs+Ik53rbAKG10cp9WSgUDRYHdtoloU1arVY7SpbMHKKSUyqeUcgdaAEv/s81S4P1HV8ArAre11pdtHdRCEjxepZQPsAh4L5me9T0twePVWufTWufVWucFFgAfJ9NiD4n797wEqKKUSqGU8gAqAMdtnNNSEnO85zE+zaCUygoUAc7YNKXtWK1WOcQZvtY6Win1KbAa44r/LK31UaVU50frp2HcuVEfOA08wDhjSJYSebwDgEzAN4/OeqN1Mh11MJHH6zASc7xa6+NKqVXAISAWmKm1jvM2P3uXyN/vUOAHpdRhjC6PL7TWyXLYZKXUr0B1ILNSKgQYCLiB9WuVDK0ghBBOwlG6dIQQQiRACr4QQjgJKfhCCOEkpOALIYSTkIIvhBBOQgq+EEI4CSn4QgjhJP4fvfMj5SJiWTcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_exact, y_exact, color='g', label='$\\sin(2\\pi x)$')\n",
    "plt.scatter(x_train, y_train, s=20, fc='none', ec='b', label='train data') \n",
    "#fc = facecolor, ec = edgecolor\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Regression(object):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LinearRegression(Regression):\n",
    "    \"\"\"Linear Regression Model\n",
    "    y = X @w\n",
    "    \"\"\"\n",
    "    def fit(self, X, y):# training\n",
    "        \"\"\"y = X @ w\n",
    "        X_t @ y = X_t @ X @ w\n",
    "        w = inv(X_t @ X) @ X_t @ y \n",
    "        \"\"\"\n",
    "        self.w = np.linalg.pinv(X) @ y\n",
    "    def predict(self, X):# testing\n",
    "        return X @ self.w\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polynomial Feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import functools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PloynomialFeature(object):\n",
    "    def __init__(self, degree=2):\n",
    "        if not isinstance(degree, int):\n",
    "            raise TypeError('degree should be int')\n",
    "        self.degree = degree\n",
    "    \n",
    "    def transform(self, x):\n",
    "        \"\"\"\n",
    "        Params\n",
    "        ======\n",
    "        x: ndarray (sample_size, n)\n",
    "        \"\"\"\n",
    "        if x.ndim == 1:\n",
    "            x = x[:,None]\n",
    "        x_t = x.transpose()\n",
    "        features = [np.ones(len(x))]\n",
    "        for d in range(1, self.degree+1):\n",
    "            for item in itertools.combinations_with_replacement(x_t, d):\n",
    "                features.append(functools.reduce(lambda x, y: x * y, item))\n",
    "        return np.array(features).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 1.]\n",
      " [1. 2. 4.]]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1,2])\n",
    "feature = PloynomialFeature(2)\n",
    "X = feature.transform(x)\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzddXhURxfA4d+NC4GEJLhLcA+uheIuRQoUijsUJ8W9/YAiLV4cirtDIWhwd4IHCQ6BeHa+P4ZCQpxsdiPzPs8+kHtn7z2blpPJ3JkzmhACRVEUJekzMXYAiqIoimGohK8oipJMqISvKIqSTKiEryiKkkyohK8oipJMmBk7gKg4OTmJbNmyGTsMRVGUROPs2bMvhRDOEZ1L0Ak/W7ZsnDlzxthhKIqiJBqapj2I7Jwa0lEURUkmVMJXFEVJJlTCVxRFSSYS9Bi+oihJR1BQEF5eXvj7+xs7lCTBysqKTJkyYW5uHuP3qISvKIpBeHl5YWdnR7Zs2dA0zdjhJGpCCF69eoWXlxfZs2eP8fvUkE4icu0aHDwIz58bOxJFiT1/f38cHR1VstcDTdNwdHSM9W9LqoefCFy9Cj//DN7ekC0bXLoETZvCrFlgbW3s6BQl5lSy159v+V6qHn4C9/IlVK8OXbvC3btw6BDcuwc+PtCpk7GjUxQlMdFLwtc0bZGmac81TbsSyXlN07SZmqZ5app2SdO04vq4b3KwaBHUrAkdO4KpqTxmbw9LlsDevfKHgKIoSkzoq4e/BKgVxfnaQO5Pry7AHD3dN8k7fhzq1Qt/3NoavvsOTp40fEyKoiROehnDF0Ic1jQtWxRNGgLLhNxe64SmafaapqUXQjzVx/2TMnt7ePYs4nPPnkGqVLG/5vOPzznhdYIbL29w4+UNvN574R/sj3+wP5qmkS5FOtLapiVrqqyUzlSaUhlLkdIyZdw+iKIkYOXKleP48eORnn/06BE//fQTz549w8TEhC5dutC3b18DRqgfhnpomxF4FOprr0/HwiV8TdO6IH8LIEuWLAYJLiFr3Rr69oX27cHW9svxEyfg1i34/vuYXeeS9yXWXl3LLs9dnHt67vPxdCnSkTVVVmzMbXCycSJEhPDg7QNOep3E+6M3ABoaRdIVoXn+5rQo2IIcDjn0+AkVxfiiSvYAZmZmTJ06leLFi+Pj40OJEiWoXr06+fPnN1CE+mGohB/R4+QIN9MVQswH5gO4urom+w13v/8eypWDChVg0CDIkQP+/RdmzoSFC8HCIvL3+gX5sfbqWuaencsJrxOYaqaUzVyW8d+N57vs35HfOT/2VvaRvv+t/1tOPz6Nh5cHuz1343bADbcDbpTNVJYBZQfQKG8jTE1M4+FTK0r8+PjxI82bN8fLy4uQkBBGjBhBixYtSJEiBVeuXKF27dpUqFCB48ePkzFjRrZs2YK1tTXp06cnffr0ANjZ2ZEvXz4eP36MhYUFZcuWxdbWFnt7ex4+fIiDgwPnz5+nYcOGuLm5Ub16dYYPH8779++ZOXOmUT+/oRK+F5A51NeZgCcGuneipmmwYAFs2CAf1D5/DkWLyqRfsGDE7wkKCWLxhcWMPTSWxz6PyeuUlz9q/kHbwm1xtHGM8b3treypnrM61XNWZ2TlkTx895A1V9Yw/9x8mq1rhoujC0PKD6FdkXYq8Sux0m93Py48u6DXaxZNV5TptaZH2Wb37t1kyJCBHTt2APDu3bsw52/fvs0///zDggULaN68ORs2bKBNmzZh2ty/f5/z589TunRpUqZMSYUKFejfvz8VK1akSpUqzJo1i5QpUzJmzBhGjhzJ8+fPOX/+PFu3btXr5/0WhpqWuRX46dNsnTLAOzV+H3OaBs2awfbtcOoUzJ8febLfenMrBWYXoOv2rmRJlYX9bfdzrcc1+pXpF6tkH5EsqbIwqPwgbvS8wdpma7GzsKPj1o64LnDl6MOjcbq2ohhCoUKF2L9/P0OGDOHIkSOk+uohWPbs2SlatCgAJUqU4P79+2HOf/jwgaZNmzJ9+nRSppTPta5evUrBT/8gb9y4QZ48eQCoVKkSQgimTZvG6tWrMTU1fqdILz18TdP+AaoATpqmeQGjAHMAIcRcYCdQB/AEfIGf9XFf5QvvD9703tWbddfWUcC5ANtabaNu7rrxstDF1MSUHwr8QLP8zVh3bR0D9g6g4uKKtCnchhm1ZpDaOrXe76kkLdH1xOOLi4sLZ8+eZefOnQwbNowaNWowcuTIz+ctLS0//93U1BQ/P7/PXwcFBdG0aVNat25NkyZNAPDz88Pf3x8HBwcePXqEo6MjFp/GWS9fvszTp09xcnLCzs7OQJ8wavqapdMqmvMC6KmPeynhrbmyhu47uuMb5MuEqhMYVG4Q5qYxL6j0rTRNo3mB5tTNXZfJRycz+dhkDt47yLLGy6iavWq8319RYuvJkyekTp2aNm3akCJFCpYsWRKj9wkh6NixI/ny5aN///6fj1+7do18+fIBcP369c9/f/r0Ka1bt2bLli306dOHPXv2ULNmTb1/nthSK20TsYDgAHrt7EXLDS3J45SHC90u4FbRzSDJPjRbC1vGVR3HiY4nSGGRgu+Xfc/gfYMJ1gUbNA5Fic7ly5cpVaoURYsWZcKECQwfPjxG7zt27BjLly/nwIEDFC1alKJFi7Jz584wwznW1tacO3eOa9eu0aRJE6ZOnUq+fPkYMWIEo0ePjsdPFXOa7HwnTK6urkJtcRixB28f8MO6Hzj95DQDyg5gUrVJBk/0EfEN8mXAngHMPTuXatmrsfaHtWqIRwHC9oAV/Yjoe6pp2lkhhGtE7VUPPxE6++QspReW5uarm2xsvpEpNaYkiGQPYGNuw5x6c1jccDFHHh6h9MLSXH9x3dhhKYqCSviJzs7bO6m0pBJWZlac6HiCxvkax+v9QkJgyxb48Udo2BCmToU3b6J/X/ui7TnY7iA+AT6U+bsMRx4cidc4FUWJnkr4icii84to8E8D8jrl5USnE+Rzjt9fj4ODoXlzGDNG1u356Sc4dw6KFIE7d6J/f7nM5TjV+RQZ7DJQY0UNdt7eGa/xKooSNZXwE4n5Z+fTcWtHquWoxqH2h0iXIl2833PRIlme+cQJ6NxZ1uBfuRL69YNu3WJ2jSypsnC4/WHyO+en4eqGrL6yOn6DVhQlUirhJwJzTs+h6/au1M1dl60tt5LCIoVB7rtkCQwdGr58Q8+ecP48eHnF7DrOts4c+OkA5TKXo/XG1irpK4qRqISfwM09M5ceO3tQ36U+G5pvwNLMMvo3xZAQsuLm69cRn3/5EiLaLtPSEjJkgFevYn6vVFap2NV6FxWyVKDNxjZsubHl24JWFOWbqYSfgK29upYeO3pQz6Ue65uv12uy37BBlmcoUEBum1itmhyfD+2/mj1f8/KSr5w5Y3dPG3MbtrfajmsGV5qvb84ezz3fHL+ifKuZM2eSL18+HBwcmDx5MgCbN2/m2rVrn9ssWbKEJ0++lPvq1KlTmPOJlhAiwb5KlCghkqt9d/YJ87HmouKiisI30Fev1167VohMmYTYt08InU4If38hFi4UwtlZiGvXvrQ7flyIdOmEOH/+y7F374SoWVOIQYO+/f6vfV+LonOLCuvx1uLEoxPffiElUbkW+n8uI8qTJ4+4e/dumGPt2rUT69at+/x15cqVxenTpw0dWqxF9D0FzohIcqrRk3pUr+Sa8E8/Pi1sJ9iKQrMLiTd+b/R6bZ1OCBcXIQ4eDH9uwgQh2rcPe+yff4RwchKicmUhGjUSwsFBiK5dhQgMjFsc3h+8RY4ZOYTz787C85Vn3C6mJAoJIeF37dpVmJubi4IFC4pp06aJnj17imPHjgkHBweRLVs2UaRIETF58mRha2srXFxcRJEiRYSvr2+YHwC2trbCzc1NFC5cWJQuXVo8e/ZMCCGEp6enKF26tHB1dRUjRowQtra28f55YpvwDVUeWYmhR+8eUW9VPZxtndnTZk+U9eq/6fqP4P17qFw5/LlWrWTd/dBatoRGjeTQzsePMGsWZMoU9zjS2KZhV+tdlP27LLVX1sajo0ecq3kqiUi/fnBBv+WRKVoUpkddlG3u3Lns3r2bgwcPsn37dkDudtWgQQPq1atHs2bNANi1axdTpkzB1TX8gtWPHz9SpkwZJkyYwODBg1mwYAHDhw+nb9++9O3bl1atWjF37lz9fjY9UWP4CcjHwI80WN0Av2A/dvy4g/R26fV+D3NzCAyUC6q+5ucnz3/Nygrq1pVz8vWR7P/j4ujClpZbePjuIQ1XNyQgOEB/F1eUeGJhYUG9TxtNhy6h7OHhwQ8//ADAjz/+aKzwoqR6+AmETuhou6ktl7wvsb3VdvI7x8/WaenTg4sLrF8ve++hzZsn59obUoUsFVjaaCktN7Sk967ezK8/37ABKMYRTU88ITM3N/9cdtzU1JTg4MRTJFAl/ARi1MFRbLqxiek1p1M7d+3YXyAgQNY8ePMG3r6V4zY+PvDhA/j6yu67vz8EBbE2dxAbO+g4tciEfAVMCNSZcfy8DfZXrek/MgVscwQnJ0iTRnbpLfU3OygiLQq24MKzC0w+NpkS6UvQ1bVrvN5PUb5mZ2eHj49PpF/HRJkyZdiwYQMtWrRg9eqEudZEJfwEYPut7Yw/Mp4ORTvQp3SfsCd9fOD+fTn4/ugRPH4sJ88/fQre3vDihXx9/Bjj+2U2NaW3ZoJuv0Dbp8MOHfWB+gD9InhDunRyQn7+/HIuZ+HCULIk6HFTh/FVx3PB+wK9d/WmYJqClM9SXm/XVpTotGzZks6dOzNz5kzWr19P+/bt6datG9bW1nh4eMToGtOnT6dNmzZMnTqVunXrhttNKyFQ5ZHjkb+/HBOPameze2/uUXx+cYqZZmJ3gUlYXL8FV67AjRuyYM3z52HfoGmy550+PaRNC87O8uXoCKlTg4MD2NtDypTyZWsrXzY2sqduZiavEVpIiAzW11f+ZvD6tVx15e0tf8g8eACennD1qjwOYGIChQrJp7y1a8tiOzY2cfp+vfF7Q8kFJfkY9JHzXc8bpHyEYjhJvTyyr68v1tbWaJrG6tWr+eeff9iyJX4XGMa2PLLq4ceD9ethwgSZH83N5cPOCRPk6tTPHj4kaMc2ri8Zwdk778nx6gqf+tiyR503LzRoIFc3Zc8OWbJA5szynJme/7OZmn75weDsHPWKqufPZV0FDw84fhwWL4a//pJPdqtV+1JW09Y21mE4WDuwueVmSi0oRZuNbdjTZo/aHF1JNM6ePUuvXr0QQmBvb8+iRYuMHVJ4kc3XTAivxDgPf+FCIbJnF2LXLiFCQoR4/lyIoUOFyJFDiNce14VwcxMif365BAKElx3iSc1yQvz+uxAHDgjx4oWxP0Ls+PsLsXevEP36CZEli/xctrZCtGsnxJkz33TJhWcXCkYjxrqP1W+silElhHn4SY1aeGVEAQFCpE0rxKVLXx38+2/h6VRKfrtNTISoVk1cHNxO5O2JGLwnDktWE5qQECEOHRKiUychUqSQn7d8eSHWr5fnYkin04nWG1oLkzEmwv2eezwGrBiSSvj6F9uEr+bhx5BOJ0cuChSQQ9UFC8Ls2fL4f86elcM2hQohx8NnzpTDIx07ktbOl5lZp8LjxzzasJgqqbdiV6Qk46tNMNpn0jsTE6hUCRYskMV2/vgDnjyBZs2geHHYtg2EwN8fli+Hvn1lrX1Pz7CX0TSNOXXnkCt1LlptaMVL35fG+TyKksSohB9D3brBqlUwZ46cFDN7NqxYAd27f2mjaXKchm3bIF8+mdGyZ4fduzk5/xJrMvYnJI0zbTe1JUgXxKqmqxLM1oR6lyqVXE15+7b8Rn34AA0a4F+iHK1ynWbFClm07d07KFtW7qQVmp2lHWuareGl70u6bOsifx1VFCVuIuv6J4RXQhnSuXhRiIwZhfjwIexxHx8h0qf/MoQTeM9LbLdsIocyChQIU7CmdWshJk0SYsLhCYLRiKUXlsZbvP8VRNPp4u0WsRcYKMT8+eKFeTqh0zRZkOfVKyGEEF5espjbsWPh3/bb0d8EoxGLzy82bLyK3qkhHf1TQzrxYPNmOfnk64knKVLI+jNbtgDu7piXLEpN3U4mp5rEllHnCK5QhSdPoH9/OHUKyjW5wCj3UbQs2JK2hdvqPU4h5LCTi4uMzdkZhgyJ1RT9+GNuzqXSnfku3Q1E7z5y2CdvXti8mYwZ5fdofgSLbAeUHUClrJXos6sP997cM3zcitEIISeE7dwpZwbH1du3b5k9e/Y3vbdOnTq8ffv2m++dIkXUmxbFJbbYUAk/BoKDw+/69B9LC0HxY7Pg++/ByQmzS+fJv2wok6dZYGkpc5qvL+x3D6D3gXY42TjxV52/Pi/NvnhRjmOPHClnOcZl5GLIEFi6VI6gBAbCyZPw8KGsgxNR7RxDe/AAshRKhcmM6bL4fqZM0LgxdOhAsZzvI/xHbWpiyrJGy9A0jbab2hKiSwAfRIl3165BiRLy8c/MmeDqKv/+7t23XzOqpBoSzT+QnTt3Ym+v30KGoRkq4Rt92CaqV0IZ0jl2TIjcucOXBA70DxErUnaTQzgNGshi8aGEhHwZVnHb7yYYjdh+c/vncx07yqGMwYOF+PVXIXLlkpfx84t9jI8eydLFn0ZJwsRQsqQQmzfH/pqx9eyZEP/7nxC9ewvx559CvPmqsvONG3II7PP3MSBAfnATE/HGIZsY1+RspNdeemGpYDRi6vGp8fcBlHgV0yGdN2/kEOrChV8md/n6CtGlixC1a3/7/Vu0aCGsrKxEkSJFxMCBA8XBgwdFlSpVRKtWrUS+fPmEEEI0bNhQFC9eXOTPn1/Mmzfv83uzZs0qXrx4Ie7duyfy5s0rOnXqJPLnzy+qV68ufH3D71dx9+5dUaZMGeHq6iqGDx/+uVSyj4+PqFq1qihWrJgoWLCg2PzpH+bXsUXW7mtqWmY80OmEqFdPiCZNhPhv34Q7njqxK3t3IUDoBg2OctrhiUcnhMkYE9Fhc4fPx2bOFKJChbDPBQIDhWjcWIghQ2If44IFQrRpE/G5WbPkD5f4tHmzEKlTC9GhgxB//CFEy5ZyQ5WjR8O2q15dLkUI/Xzh/qpj4pFJZhFiYSnE4sURXl+n04l6q+oJq/FW4tbLW/H3QZR4E9OEP326EK1ahT8eFCQ7SBcufNv97927JwoUKPD564MHDwobG5swm6G8+tRj8vX1FQUKFBAvX74UQoRN+KampuL8p12BfvjhB7F8+fJw96pfv75YulQ+p/vzzz8/J/ygoCDx7lPH8MWLFyJnzpxCp9OFiy2ydl9TY/jxQNNg3To5w7JkSXBMLdhbsD+17s0huP9gtN8myymJEQgIDuDnLT+T0S4j02pO+3x89myYPDnscwFzc5gyBRYuhKCg2McYl+GguPD2hg4dYO9e+PtvOTnnn3/k1MumTWXdtv8sXy7bFS0KgwfLZyNFe5Tj+MyzmFQoDz//DD16hPsGaJrG3LpzsTS1pNO2TuiEDiVpOnUKatUKf9zMDKpXl+f1pVSpUmQPtXHzzJkzKVKkCGXKlOHRo0fcvn073HuyZ89O0aJFgbDlkUM7duwYrVq1AqBt2y/P64QQuLm5UbhwYb7//nseP36Mt7d3uPfHtF1sqYQfQ1ZW8PvvsmbZww6j6eY/Hfr0wWzK5PC1aUKZeGQi119eZ379+aSy+lJM6d49OTX9azlyyLn9sR2rrFVLPtz6ekPykBBYtkxWaYgvS5dCkyZyzDW0mjVlYt+8+cuxtGnlP9hp02TZnypVZMmg5j2dYc8eGDRIzn2tX18WjgslY8qM/FHzDw4/OMzs0wYY71SMwsFB1giMiJeXLBmlL7ahelzu7u7s378fDw8PLl68SLFixfD39w/3HstQ1WOjKo+sRZAXVq5cyYsXLzh79iwXLlwgbdq0Ed4jpu1iSyX8WDLftRXbqWOhfXtZ0zuKZH/l+RUmHZ1Em8JtqJUrbJclRw65UOtrd+7I0jaxfT6UMSN06iST7H8Pf2/f/jK7qG7d2F0vNh4+/LTYLAKFC4efYaFpsuzOsGHQpUuof8BmZvKn6sKFsH+/3Jbr6dMw721ftD01c9Zk6P6h3H97H5B136ZMgSJFZLmhxo3h6FH9fkbFcNq2hblzZaXv0M6ckbN26tT5tutGV/L43bt3ODg4YGNjw40bNzhx4sS33QgoX7785xLJK1euDHOPNGnSYG5uzsGDB3nw6R/H17FF1i6uVMKPjbt3oV072TWfMyfKZB+iC6HT1k6kskrFHzX/CHe+Rw85q+bDhy/HAgNhwACZuL+lPtpvv8kRkfbt5fvLlZPrvnbsiLpiZ1y5uMgZQRE5eRLy5InlBTt2lIvXbt2Sq7Lu3Pl8StO0z5uk9NzZk4AAQd264O4up6QePSp/2/nhBzkMpyQ+pUvLzXlKl5ZTdd3d5Uy2OnVkX8Da+tuu6+joSPny5SlYsCCDBg0Kd75WrVoEBwdTuHBhRowYQZkyZb75M8yYMYO//vqLkiVL8i7Ur+utW7fmzJkzuLq6snLlSvLmzRthbJG1i7PIBvcTwiuhPLQVQsipM8WLC2Fv/+XJbRSme0wXjEasvLQywvMhIUJ07ixEhgxCDBggxLBhssBaw4Zy0VRcBQYabuHV69dyo/N9+8IeX7lSiKxZ47Dh+enT8klwxoxC3LwZ5tS049MEoxE9Z68VlSqFf2Z+6pQQ6dLJiUBKwhCbhVc6nRB79gjRvLkQlSoJ0aOHEFevxmNwiZSapRNfevaU366tW6Nt+ujdI5FiYgpRe0XtCJ+sh3bpkhBjxwoxapQQHh4JbHVsLBw6JESaNELUqSNnGVWuLItnhikk9y0uXpQ/TdKnF+L69c+Hg0KCRPF5xYWlWzqxcMWbCN9aqpQsQKokDGqlrf6pWTrx4cQJOV7Qt698mBiNX/b8QrAumD/r/Bnhg5vQChWCESNg9GgoUybKUaIErVIl+SC6VSu570rv3rIoWmRj+zFWuLD8nV6nk094b90CwMzEjPn15hNg/pyN79wifGvKlGFnCClKcqcSfnSCg+WAe8aMMG5ctM133d7F+mvrGV5xODkcchggwITDxgbatAE3Nzkd01xfdeEKFPiS9KtXl1M1gBIZSlDOtA87n8/lpFfYhwje3nD6tPwhqiQcwlhzh5Ogb/leqoQfnTlz5NSA6dOj3cPVL8iPXrt6kccxDwPLDTRQgMlE3rywe7eculG9+uetFld2HIuJb3qaLe7Jh49yefytW3KmTvfu+p3Cp8SNlZUVr169UklfD4QQvHr1Cisrq1i9T21xGJWnT2H4cDnXsWnTaJtPOjqJu2/u8u9P/2JpZhlteyWW/qupX6uWfLm7ky2DHVNrTOWXY61IU3sB6R9348MHOfo2dKixA1ZCy5QpE15eXrx48cLYoSQJVlZWZMqUKVbvUZuYR+Wnn2DtWrmpeK5cUTa98/oO+Wfnp1n+ZqxssjLKtkocbdsGjRrJxQWbNiFMTKi2rBrnn15gV52bFMvjjKX6easkU1FtYq6XIR1N02ppmnZT0zRPTdPC9as0Tauiado7TdMufHqN1Md949XNm7LsZN++0SZ7gP57+2NhasH/qv/PAMElc/XryxKK27bBoEFomsafdf7kQ5APC+8PizLZv38P8+bJ9Q7Tp38eGVKUZCHOCV/TNFPgL6A2kB9opWla/giaHhFCFP30GhvX+8a7iRNlPYUBA6JtuttzN1tvbmVEpRFksMtggOAUevaEPn3kNopz5pDfOT/9Svfj7/N/c+pxxMVWPDwgd25ZyyddOvloxsVF/txQlGQhsvmaMX0BZYE9ob4eBgz7qk0VYHtsr220efienkKYmgrxyy/RNg0IDhAus1xE7pm5hX+QHlZMKTEXHCxE3bryv9WBA+K9/3uRbko6UXpBaRGiC7sSy89PTuXfvj3sJU6dkmu7vL0NGLeixCPieR5+RuBRqK+9Ph37WllN0y5qmrZL07QCkV1M07Qumqad0TTtjNEe7kyaJGsTRLD8+mszTszg1qtbzKg1Qz2oNTRTU1mW08UFWrTA7vlbJlWbxMnHJ1l1eVWYplu2yI3nv64pVLKkfBywbJkB41YUI9FHwo9oqdDXT4LPAVmFEEWAWcDm8G/59EYh5gshXIUQrs7OznoIL5YePJDlHzt1gvTpo2zq/cGbcYfHUc+lHrVz1zZQgEoYdnawcaOsoNa0KT/laY5rBleG7B/Ch8AvhYru35fF1SJStKhcNKYoSZ0+Er4XkDnU15mAJ6EbCCHeCyE+fPr7TsBc0zQnPdxb/37/XS53HTIk2qYjDo7AL9iPqTWmGiAwJVJ588of0qdPY9K3HzNrzeSJzxMmH538uUmuXHIhVkROn5Zj+4qS1Okj4Z8Gcmuall3TNAugJbA1dANN09Jpn2oMaJpW6tN9X+nh3vr1/r1MHK1byzq7Ubj47CJ/n/+bXiV74eLoYqAAlUg1bizrLS9YQFl3T1oXas2U41M+l1CuX18WO12zJuzb3N1h1y5ZkldRkro4J3whRDDQC9gDXAfWCiGuaprWTdO0bp+aNQOuaJp2EZgJtPz0cCFh+ecf+PgRunaNspkQgv57+2NvZc/Iygl/hmmyMXasLOrTvTtTsnXBRDNh2L/DALkJ/bZtMHCgXKg7YoQcu//hB/lDwNHRyLErigGohVdhbyiL0l+8GGUVs603t9JwdUNm1Z5Fr1K9DBefEj0vLzlYnzUrYybVYvSJSXh09KBMJllUJyAANm2SyyyyZJEJP0UKI8esKHoU1cIrlfD/c+6c3KNv5kxZ6jESQSFBFJhdAFMTUy51u4S5qb4qhCl6s3UrNGxIYK8eZM22kez22TnW4Vi0lUsVJSmI95W2ScKCBXKhVZs2UTabd3Yet1/f5n/V/6eSfULVoAH06YPFn7NZYtECDy8P1l9bb+yoFMXoVA8f5Lh9+vTRTsh+5/+OXLNyUTBNQQ78dED1GBMyf38oWRLx6hWV+zvgZe7H9Z7X1VoJJclTPfzorFkDPj5yR+0o/HbsN176vmRK9Skq2Sd0VlawfDnay5esPejMvTf3mH16trGjUhSjUgkfYNEiOZe7fLISm1oAACAASURBVPlImzx694g/TvxB60KtKZGhhAGDU75Z0aIwejTpdh5i0vNCjD8ynrf+b40dlaIYjUr4T57AsWPw449RzswZcXAEQggmVJ1gwOCUOBs8GMqWZeCq+1h7vw6zGEtRkhuV8Ddtkn9GscHJZe/LLLu4jN6lepPVPquBAlP0wswMli3DLCCIrUcyMePEdB69exT9+xQlCVIJf8MGOZyTP6KKztKwf4eRyioVwyoOM2Bgit7kygXjxlH8tBdNLwUzyn2UsSNSFKNI3gn/5Us4dAjRuAnnz8OhQ/D2qyHeQ/cPseP2DoZVGEZqa7VBaqLVrx+ULMncvRZsO7aYK8+vGDsiRTG45J3wt2wBnY5m/zSlRQv49VfInl0O+4aEyBIKQ/YPIaNdRnqXinwxlpIImJnBokXY+gbz115zfj3wq7EjUhSDS9YJ/+OyDTwwycZPfxTj5k04ehRu3JDVE3/9FTbd2MTJxycZU2UM1ubWxg5XiauCBdF+/ZXmF4MI3L6VYw+PGTsiRTGo5Lvw6u1bgh3TcLJUH8p7TAlz6tkzyJs/mLSjC2FqqnGp+yXMTMziJw4lUu/eyQXQ/21BWK+eXCqRKlUcLhoQgK5IYR69uEPHCaXY11WVXFCSFrXwKiI7dmCmC8KhY5Nwp9KlA8eqy7n15gbjq45Xyd4Inj+HUqXg7Flwc5O/cZ0/L495e8fhwpaWmMybT9bXIVRd6cEuz116i1lRErrkm/A3buSlRQY8ncqEO+UfFMCD7KMp6FCSxnkbGyE4ZcQIqFNHVqyuWRNq1IBVq+QWhSPjWpG6cmV0P7Vl0HFYuPwXdEKnl5gVJaFLnkM6QUHg6MiN4q3oxjz275fP9P7z87wZLHnWj71t9lE95/f6v78SJSEgZUq4fVv+thXa06dyC9v376NcJxe9Fy8IyJ2Dk6k+4LV1OT8WibponqLEl9u3YedO+f99nTry/++4UEM6Xzt9Gnx8yNWtOra2UK2aXH91/Dj0H+bD0rsTKJG6qkr2RhISAn5+kCZN+HNp08pzwcFxvImzM+ZTplHpIZyfMoCgkKA4XlBRYkeng549ZUWX69flHg0VK8r9l0JC4ueeyTPh798PmoZZ9e/YvBnat4fZs+GXX+CEmIGwecFfjScaO8pky8xM7mGyZ0/4c3v2yHPmeqhMbdKhI2+KuDBg03NWHp0T9wsqSizMmgUXLoCnJ8xNNYQ5NgO44ym4fh3++COebiqESLCvEiVKiHhRsaIQEVz7le8rkXJSStHwn4bxc18lxtatEyJ7diGuXPly7OpVeWzdOv3dR3fqlAjREPMqpxB+QX76u7CiRCNnTiFOnhRCPHgghIWFEJ06CSGEOHtWiCxZhNDpvu26wBkRSU5Nfj38Dx/AwwO+Dz9c879j/8MnwIdx340zQmBKaM2ayT3Jq1aF0qXl67vvYOhQeU5ftJIledqiLh0Of2DtWlVyQTGMkBC4exdKlgTGj5cHR4wAoHhx+awqMFD/901+8w0PH5YDwNWrhzn87MMzZp6aScuCLSmUtpCRglNC69wZfvoJTp6UD7TKlAHLeNi/JOOsJbzflp5so/7gQ9PhpLC00/9NFCUUU1M5IcFztye5Fy2CHj3kJsvIsXx7e7Cw0P99k18Pf98+uTnGV7XvJx+dTEBwAKOrjDZOXEqELC2hUiWoXDl+kj0ATk68GtaPSp5B7Pm9azzdRFHC6tIFHncejbCwkItNkH3RYcNkZyc+1gMmv2mZhQrJH6379n0+9OjdI3LNykWbQm34u+Hf+r2fkjgEB3M/R2r4+BH7u0+wT5U2xm8VQk6tCwiQhVf18UBZSfoCz1/FrHghFjsO4vmA3zAxgRUrIFMmOWvQyurbrqumZf7n2TO4ciXc+P24w+MQQjCyclxX9CiJlpkZwf/7nWyvdZwa3DrGbzt4UM4aqlYNfvgBsmWDefPiL0wl6bAYPxItpR0uCwfz/Lkct58+Xc7J/9ZkH63InuYmhJfeZ+msWCEECHHmzOdDnq88hdlYM9FzR0/93ktJlE6XSC/eWSJe3r0abdtz54RwdhZi82YhQkLksQsXhHBxEWLhwngOVEnczp2TuWjUKL1fGjVL55P9+yF1aihW7POhsYfHYmZihltFNyMGpiQU9n8uxDoIPHu2irbtb7/JGj8NG4LJp39JRYrIX8vHjYu/xTNKEjByJDg4yMU/BpS8Er67u5zb9+lf542XN1hxaQU9XHuQwS6DcWNTEoRcZeqwv25eSu6+xIvj+6Nse+gQNAlfe4+SJeXDt4cP4ylIJcELDJR1oPr3lz/8794NdfLkSdi+HQYOjGPp19hLPgn/2TO4fx/Klft8aMyhMVibWTOkwhDjxaUkOHmmr+CtFbzq3THKdra28Pp1+OOBgfDxozyvJD/370PBgrBwIaRPD69eySqvU6d+ajByJDg5QZ8+Bo8t+ST8Eyfkn2Vkdcwrz6+w5soaepfqTRrbCIq2KMlWjhwl2NuqFHnPPeTFppWRtmvZEmbMCH982TI5ahhRLSAl6WvTRk65/PdfGDRIPoi9cEH+v3J5zlHYuxeGDIEUKQweW/KZljl0KEybJsssWlnRbG0z9t7Zy72+93C0cdTPPZQk4+7Tm1AgL8Hmjmzs402zFqbkyhW2zZs3sthV0aJy3UyKFLB2rdy0Zc8eeVxJXq5ckRUv792Ti6tCmz4davxWlfxchzt3wMYmXmJIttMyfX3ltLnDh0F33EN2u6ysuPDsAhuub6BfmX4q2SvhvH0LbZvlYULBarg8f4Wz+x+ULRu+Dr+Dg9wWM29emfCbN5c/BI4fV8k+uXrwAAoUCJ/sASoLd/I/Oyg7n/GU7KOTJHv4QshqcxMnQp48EBIQzMGzKXlQswt5d0+n0epGuN93536/+9hb2cdD5Epi1qaNrMc/dMJDnhbOhouvNSHnX1Cxpg2TJkGjRsaOUEmobt2Sq8IfPPiqNIIQPMxZhZTet7F/eQes42+P7GTXw1+0SD4wOXkSjh2DE/MvYY0fs06X4a+NZ9lycwsDyg5QyV4J5+VL2LEDJk2CLA5Z8OjdGIfXvpgsHMWoUTBHVVFWouDiIhfzjxkjO57/8Vp2gCz3DvO+p1u8JvvoJLkevk4nv+nLl0PZsp8Ozp4NPXuybdZ9fr7VE13G49zvd5+Ulin1H7SSqF24AO3awcWL8usnPk84VyozVR+a8eTAE2q1dsTT07gxKgmbt7fcijMkBGrVAq9Hgp6rK1Io5X1sn3jG4zJaKVn18F+/ltOgyoTeqtbDA9Klw7b6U1457mBQuUEq2SsRypxZzp9/905+ncEuA1f7/oiVXyDvRg8hZ07jxqckfGnTwqlTcmGerS00d9hHmZBj2I53i/dkH50k18P39ZXT4R4/DrWmIXduKFSIijV9OX7vLG/H3MVOlcBVItGmDdjZwV9/yTV6zz48Y1+lTPxwUXB43j1qdMpi7BCVxEIIufbn8WNZYS/eSr5+kax6+DY2UK9eqC3CXrwAT0/u503P0Wd7qGg6SCV7JUp//gmXL8uZNiNHwoRh6Rjn1AE0HaX/7Wfs8JTEZO9euQbIzc0gyT46Sa6HD/Do0Zca6n1zbKPYqAbUb16C3Tkecb//XTI6qyWQStSEkAtnDh2Sv5ZXrf8cj24Z6X08BJNLl+XcO0WJihF695DMevggx2HPnJFTMq8sPEGwZsq/uc4ytuZgleyVGNE0WUV73Dg5bbpUgTS87dcdH3PB+0F9jR2ekhgksN496Cnha5pWS9O0m5qmeWqaNjSC85qmaTM/nb+kaVpxfdw3Ko6OcueYtrk8uJPVhpQOaelbvnt831ZJwnrWHsnMihak3PWvfCqnKJERQs7NzJwZfv7Z2NF8FueEr2maKfAXUBvID7TSNC3/V81qA7k/vboABpvN/DBPOhbk9mFohaHYmBtndZuSNDjZOBHSpxfPbeDDQNXLV6Kwb5+cHRjD3v37gPcYYnhdHz38UoCnEOKuECIQWA00/KpNQ2DZp/r8JwB7TdPS6+HeURJC0Lb0Y1bVTE/XEmqvUiXueldzY9p3lqQ4ckIO8ivK14SA0aNj1btvt7kdNVbUiPekr4+EnxF4FOprr0/HYtsGAE3TumiadkbTtDMvXryIU2ACQcsCLZn8/WSszY23uk1JOhxtHLHu9QsPU4LvoH5hl1MqCsiNljw85JhyDHr3556eY/ONzVTMUhEtPnYuD0UfCT+iCL/+VxCTNvKgEPOFEK5CCFdnZ+c4BWaimdC9ZHd+KvJTnK6jKKH1qTyY37+3wub8Fdi61djhKAnJf2P3mTJBhw4xesto99HYW9nTt3T8DxPqI+F7AZlDfZ0JePINbRQlUXCwdiBNj0HcSg1+wwbJeh6KAnDggCzgNXRojHr3px+fZtutbQwsO5BUVvG/+5U+Ev5pILemadk1TbMAWgJfd3u2Aj99mq1TBngnhHiqh3srilH0Kd+f32vYYH39NqxbZ+xwlIRi7FjIkAE6Rr1b2n9GHxpNauvU9C7dO54Dk+Kc8IUQwUAvYA9wHVgrhLiqaVo3TdO6fWq2E7gLeAILgB5xva+iGJO9lT1ZuwzmchrwdxssN7FVkjd3d7n5xpAhMaqZc8LrBDtv7zRoba8kudJWUQzhfcB7enTOwIrlH2HxYmjf3tghKcZUtSpcvy53LI9BCeRaK2px9ulZ7vW9RwoL/W13mOxW2iqKIaS0TEmBTm6cSQ8BI93k7uVK8nTkiNxeb/DgGCX7Yw+PsefOHgaVG6TXZB8dlfAVJQ56le7N77XtsHz0FP7+29jhKMYydqws09s1Zut9RrmPIo1tGnqW7BnPgYWlEr6ixIGdpR2uP//K0cwQMHYU+PsbOyTF0I4fl3PvBw2K0V61h+4f4t97/zK0/FBsLQxb20slfEWJo56lejGtjj2Wz17AvHnGDkcxtHHjwMkJunWLtqkQgpHuI0mXIh3dXKNvr28q4StKHNla2FKh3QgOZIPA8WPlLjxK8nDqFOzeDQMGQIrox+IP3j/I4QeHcavgZpTV/yrhK4oedHftzvS6qbF4+Rrx55/GDkcxlHHjIHVq6Bn9WLwQgpEHR5LRLiOdS3Q2QHDhqYSvKHpgbW5N9baj2Z0TgiZPAB8fY4ekxLdz52D7dujXT+6JGY09d/Zw7NExfq34K1ZmxtnbViV8RdGTziU682f9NFi8eY+YMcPY4Sjxbfx4uXF2nz7RNhVCMOLgCLKmykrH4jFbhRsfVMJXFD2xMrOiXusxbHWB4P/9Bu/eGTskJb5cugSbNsnefaroa+Bsu7WNM0/OMKLSCCxMLQwQYMRUwlcUPepQrAPzGqTH/P0HxLRpxg5HiS/jx8thnL7RV7jUCR0jD44kp0NOo1fuVQlfUfTIwtSCZj9OYEM+CJ42BV6/NnZIir5duwbr18uhHAeHaJtvvL6Ri94XGVV5FOam5gYIMHIq4SuKnrUt0pbFDbNg/sEXMXWqscNR9G3CBLnA6pdfom0aogthlPso8jrl5cdCPxoguKiphK8oemZmYkbrHyezpgAET58GL18aOyRFT14eu4nun9Xsy9OTWascefMm6vb/XPmHay+uMabKGExNTA0TZBRUwleUeNCiYAtWNsmFqZ8/ut9/M3Y4ih5s2QL/Vp1AkKkV9xoPwMMDXFxk3bSIBIUEMcp9FEXSFqFZ/maGDTYSKuErSjww0Uxo3+o3VhWEkFkzwdvb2CEpceDtDePaedI8eCWWfbvTZXgaVq2CVaugWTPw8wv/niUXlnD3zV3GVx2PiZYwUm3CiEJRkqDGeRuz8YcCmAQEEjxpIp6eMHw4/PwzTJmiRnoSk+XLYXqaiWgWFjBw4Ofj1auDqyts3Bi2vX+wP2MPj6VMpjLUzV3XwNFGTiV8RYknmqbRtdVUlhUB3ew5NCr1hMBAqFABrl6F/Pnl9qdKwvfx8l3K3Vkmyx+nSxfmXKFC8OBB2Pbzz87H670X478bj6ZpBow0airhK0o8qpGzBlsau6KFBHGozjh+/11ud7p4MSxdCs2bq31TEoNmtycRjJnc4OQrJ09Cnjxfvv4Q+IEJRybwXbbvqJajmgGjjJ5K+IoSjzRNwy5gGouLQqq1C+Dhw8/nateGXLlg1y4jBqhE7/598p9ewnKLTuy9kiHMqVWr5I6G9et/OTbjxAyef3zOxGoTDRxo9MyMHYCiJHUmXhVZX6MK7S66EzJ2FJYLF38+ly8fPHlixOCU6E2ciGZiQqGVQ2nwExQrBoULy579/fuyfprFp2oJr/1e87/j/6NBngaUyVTGqGFHRPXwFSWe5c8PNm+nsqA4mC1ZJruEgBBw9KhM+koC9eCBHH/r1IlSTTJx7x60bQv29nKh7e3bcgz/P78d/Y33Ae+ZUHWC8WKOgurhK0o8a98ephQsjmm/+nQ8vw1GumG1fDVTp4KZGVSubOwIlUhNmgSaBkOHAnJ/8h8jWTD7xOcJs07NonXh1hRMU9CAQcac6uErSjxLmxY2bIDDy6cy11XDYtVaGua9yYoVcjFPAprEoYT28CEsWgSdOkHmzNE2H394PEG6IMZUGWOA4L6NSviKYgCVKoHXxdx41GuLv5lgRvpBnD8PWbMaOzIlUpMnyz8/9e6jcvvVbRacW0CX4l3I4ZAjngP7dirhK4qBWFrC9D6TmF3WjKyHt6FdvWLskJTIPHwICxdChw6QJUu0zYcfHI6lqSUjKo8wQHDfTiV8RTGgDHYZ8OvbEx8LeDs4+lrqipFM/DSl0s0t2qZnnpxh7dW1DCg7gHQp0kXb3phUwlcUA+tdexSzK1phv+sAnDlj7HCUrz148GXsPprevRCCIfuH4GTjxIByAwwU4LdTCV9RDMzB2oEUg4fz0hpe9e9u7HCUr02YIJ+kx6B3v+/uPg7cO8CISiNIaZnSAMHFjUr4imIEnb8bwLzq9jgeOYPukLuxw1H+c++enHffuTNkyhRlU53QMWT/ELLbZ6dria4GCjBuVMJXFCOwMrMiu9sUnqSAl790k6uwFOObMAFMTWHYsGibrri0ggvPLjCh6gQszSwNEFzcqYSvKEbSstTPLKmXkTTnbxK4c7uxw1E8PWHJEujSBTJmjLKpX5Afww8MxzWDKy0KtjBMfHqgEr6iGImJZoLriLnctYc3/buDTmfskJK3sWNlUZwY9O5nnJzBo/ePmFJ9SoLZ3CQmEk+kipIE1chfj3UtC5L21mN8ViwydjjJ17VrsGIF9OoF6dNH2fTFxxdMOjqJ+i71qZwtcdXFUAlfUYyszojlXE4Dfm6DICjI2OEkT6NHg61thPXuvzbu8Dg+Bn7kt+8T317FKuEripEVylCUQ11qkObxW7xnTTZ2OMnPxYuwbh388gs4OUXZ9MbLG8w+PZvOxTuTzznxlTlVCV9REoBmg5dwMosJZuMnRrwjthJvQoaPJMDGnnoH+lOqlNyy9ustC/8zcO9AbC1sGfvdWMMGqScq4StKApDOLj23BvyM4xt/7o7pZ+xwko0Adw9Mt29lZYZBdBtqz/Tpcs1VqVJw/nzYtns897Dj9g5GVBqBs62zcQKOI00k4Pm/rq6u4oxaeq4kE35BfngUSU3J+4HYPHiKqXMaY4eUtAmBV+7vsHl0A/uXdzCxs/18aulSmDMHTpyQXwfrgikytwgBwQFc7XE1Qc+71zTtrBDCNaJzcerha5qWWtO0fZqm3f70p0Mk7e5rmnZZ07QLmqapDK4oEbA2tyZowjhs/XRc6d/G2OEkfXv3kunOIV50GR4m2QO0aSOHdT5tTsaCswu49uIaU2pMSdDJPjpxHdIZCvwrhMgN/Pvp68h8J4QoGtlPHkVRoEajAeypkJa8/+zj3c1Lxg4nQXJ3h0aN5AbwlSrJ3nislzDodODmhpd5NoI7dAl32tQUnJ3h3Tu5T+2IgyOokq0KDfM01MtnMJa4JvyGwNJPf18KNIrj9RQlWdM0jSzTl6DT4FbPlsYOJ8FZsEDuKVu/PuzcCYMGwezZsmx9rEan16+Hc+fYWWYs2/ZYhDt9757cXD5vXhh5cCRv/N8wo9YMtMS+PZkQ4ptfwNuvvn4TSbt7wDngLNAlmmt2Ac4AZ7JkySIUJTna3aSICNEQd/7dYOxQEow3b4Swtxfi9u2wxz9+FMLFRQh39xheKDBQvqFAAXHpfLBwdhZi1y4hdDp5+tEjIcqUEWLsWCEuPL0gTMaYiF47eun1s8Qn4IyIJL9Gu4m5pmn7gYiq+v8ai58r5YUQTzRNSwPs0zTthhDicCQ/gOYD80E+tI3FPRQlySgxcz1vd7vwtndnxOVGaCZqQt2OHXLD91y5wh63sZHFLVevjuGG8AsXwq1bsHUrhYqasno1dO8uZ+ekTg03bkDv3uDmJvhuWW9SW6dOtNMwvxZtwhdCfB/ZOU3TvDVNSy+EeKppWnrgeSTXePLpz+eapm0CSgERJnxFUcApYy4Od29KpanrOTL3Vyr2mGTskIzu40eZkCPi4CDPR8vHR66qrVwZ6tUDoGpVmeTPn4cPH6BYMbCzg38ur+bIwyPMrzcfB+sI56MkOnHtNmwF2n36eztgy9cNNE2z1TTN7r+/AzUAtZmnokSj/ITl3E9rSboxU/jw8Y2xwzG6ypXluH1AQPhzmzbFsHc/ZQo8fw6//y679J9oGhQvLh8C29mBT4APA/cNpHj64nQo1kF/H8LI4prwJwPVNU27DVT/9DWapmXQNG3npzZpgaOapl0ETgE7hBC743hfRUnyTC2t8J8wltzPg9k/LPGU4I0vefLAd9/Bjz+Ct7c85usri1zeugWtWkVzgadPZcJv0UKurIrCyIMjeerzlDl152BqYqqfD5AAqIVXipKQCcH1whlwuvuM1xdPkCdXaWNHZFT+/nJmzooVkDkzPH4M5crJRVLRbFAFXbvK3ayuX4ecOSNtduHZBUrML0GX4l2YU2+Ofj+AAUS18EolfEVJ4F4f/xf7Ct+zrnpGftj9MFHVX48v79/LqZNp0kRbzVi6fBmKFpXlj2fMiLSZTugov6g8d9/c5UbPG4ly7D7eVtoqihL/Uperxu2GFWmy/zEbN040djgJQsqUUKRIDJO9ELISZqpUMGpUlE0XnlvICa8TTKk+JVEm++iohK8oiYDL3HX4W5ri4DaG5x+8jR1O4rJ9O/z7L4wZE/k0H+DZh2cM2T+EKtmq0KZw0ixtoRK+oiQwGzfKcWlra8iaVXZK/ezS4us2iGq3glk2obmxQ0w8AgNhwAC5ZLZbtyib9t7VG78gP+bVm5f4V9RGQiV8RUlAZsyQmy4NHQovX8rO6bVrUKcOpO4/ludZnWg07zB7rm41dqiJwssxf8Ht2+yvM42X78wjbbf5xmbWX1vPqMqjcHF0MWCEhqUSvqIkEO/eyTVB+/ZBgwZyx71ChWDNGggJgQ1bzbGfu4Rcb+DioJ94H/De2CEnWEFB8MuP3phPGsOlDDWZ+6A2uXLB1Knh2771f0uPHT0okrYIA8sNNHywBqQSvqIkEP/+C2XKQPbsYY+bmEDHjnJxkUWturyqV5Xe+97x+/KohyiSs9GjofahIaQ086XwgRmsXy8n6sydC1u+Wh46ZN8QvD96s7DBQsxNI/8tIClQCV9REoiQELAIX7gRAEtLCA6Wf3ecuwzN3Jyyv//Dgbv/Gi7ARCIgAM7/eYwaT5aiDRwoV2wh5+1PnAjTp39pu/fOXuafm0//Mv1xzZD0K7erhK8oCUTlynDoELx4Efa4EHKhUe3anw5kzIjJ2PHUvQ1rx7fiQ+AHg8eakD19FMz//HrKDP9r2BqPlSrB1avy72/939Jxa0fyOeVjXNVxRojU8FTCV5QEIk0aWbWxbt0v+6m+eAH9+8vdl3788Utbi379+ZA/F79ueMGvm/sYJ+AEKt3muRQIusiHcX/IByGh3Lz5Ze5+v939eOrzlKWNlmJlZmWESA1PJXxFSUAmTpSJvWFDOWU8Z075MNfdXZYB/szMjBQLl5HRB3JMW8z2W9uNFXLC8uQJVuN+5XK66oy80CTMpigBAXIqfseOsPXmVpZeXMqwCsMombGk8eI1MFVaQVESIJ0OXr+GFCnAKorOZ0jPHmhz5lC/hwOLf7tBGttkvvF506awcyev3C9TtUsunJxkrbSPH+Hvv+V0/BkLvSm+sDDpU6TnVOdTWJhG8uAkkVKlFRQlkTExASenqJM9gOnk3wjOmJ6pa97QY0MHEnIHLt5t2iRXrY0ejWPpXJw6Jbc+9PCQ9e7/+ANWr9HRaUd73ge8Z0WTFUku2UdH9fAVJbHbvRtq12ZcJXCeOodurslwuua7d5Avn3wQcvo0mEc8vXL6ien8sucX/qrzFz1K9jBwkIaheviKkpTVqoX4qS3DjmksWtyHS96Xwpz28ZFj1y4ukC4dNG4se71JytChskj+woWRJvsLzy4wZP8QGuRpQHfX7gYOMGFQCV9RkgBt2h9oTs4s3Shos6rZ56maHz9CtWpySGP1ajh7FmrVkkl/2zYjB60v+/fLFVX9+oFrxHPpfQJ8aLm+JU42Tvzd4O8kWysnOirhK0pS4OiI6eIl5HsWTLv1t+m+oztCCBYtktMQV62SW/hlzCj3AVm7Fvr0kYu9ErW3b+Hnn+XT2PHjI2wihKDTtk7cfn2blU1W4mTjZOAgEw6V8BUlqahdG7p3p78HPNqyggXnFrB+vZzb/3WHtlIlOUX93DnjhKo3ffrIrQuXLZPlRSMw69Qs1l5dy8SqE6mSrYph40tgVMJXlKTkf/+DXLlZs90Kt409eWHlQYoUETe1tZVbBiZaGzfC8uVyNW3JiOfSH390nAF7B9AwT0MGlx9s4AATHpXwFSUpsbVFW76cNG+DWLrTkkdlmvD3uifhmt29C3fuQIkSRohRHx4/lmNTxYvD8OERN3n/mB/W/UDWVFlZ0mhJsh23D00lfEVJakqXRps4Gz6sTQAADoVJREFUkboXPtLp3GtWBjZj6vQA/Pzk6YsXoVEjGDLkq9W7iUVwMLRqBX5+sHJlhLNyfIN8abi6Ie8D3rOpxSbsreyNEGjCoxK+oiRFAwdC3bpM2aOjkM6DaZ6dSZtOkDkz1KsHXbrIJonS6NFw5AjMmSMf1n5FJ3S029yOc0/PsarJKgqlLWT4GBMoM2MHoChKPDAxgaVLMS1WjP3bPpD9p+V0XZWdngXGkCkTmCXWf/n79smCQz//DG3bRthkjPsY1l9bz+/f/079PPUNHGDCpnr4ipJUOTrC2rXYv/DhyN6M/HFyLAffLE68yf7hQ2jdWq6onTUrwibzz85n7OGxtC/aPsnvXvUtVMJXlKSsTBm02bMpdP4xa05locv2Luz23G3sqGLv40e572NAAGzYEK7sMcDG6xvpvqM7tXPVZn69+eohbQRUwleUpK5zZ+jbl6Z7HjL8dgYar2mM+313Y0cVczodtGsn9yhcvTrCcXv3++602tCKUhlLse6HdUl+q8JvpRK+oiQHU6ZAjRqMXP2Uli/SUm9VPTweJZKCOmPHyl79//4XatuvL449PEb9f+qTK3Uudvy4A1uL8L1/RVIJX1GSAzMzWLMGLXdu/l70iupvHKi9sjanH58O00wIOHVKdqRPngSjF9OdP19Wfvv5Z/jll3Cnjzw4Qs0VNclgl4F9bfeR2jq1EYJMPFTCV5Tkwt4e9u7FxNGJ9Yt9KfHelqrLqn4e3rl7Vy5YbdNGlpZv21bWIrtzJ/yl9u2DqlVlvf40aWTdspcv9Rzv2rXQrRvUqQPz5oWrD3Ho/iFqr6xN5lSZcW/nTga7DHoOIOlRCV9RkpOMGWHfPkzNzNmzTFA6KC21VtRi07Vt1Kolk/2NG7BmjfyzXTuoWRMCA79cYu1aaN9ezuV/9Ur+JhAUJOvzvHunpzj37JHBlC8P69aFW1y16fomaq2sRVb7rBxsd5D0dun1dOOkTSV8RUlucuWCPXsw++jHnvl+1NPlotm6xlDsb/r1k1P4Qf7Zpw9kygSbN8tjwcEwYIAsY9OypZwskz07/PUXFCkiO+JxtmuXrN+cP7+s4fzVcuBZJ2fRdG1TiqQtgns7d9KlSKeHmyYPKuErSnJUpAi4u2MaFMzaP72p/diV2/k70W93P4J1wWGafv89XLgg/37+vBwZKl06/CU7dfryg+GbrVkjp1/mzSvHjey/lEQI1gUzYM8A+uzuQ4M8DTjQ7gDOts5xvGHyohK+oiRXRYrAkSOYWFmzYfl1fr7clBknZ1B7ZW1e+b763OzePXD+lFdDQiLdUApz8zjW158/X9bIKVsWDh78clPgqc9Tvl/2PdNOTKN3qd5saL4BG/PEWAjIuFTCV5TkzMUFjh7FNHNG5m/YzF6vFhy6507huYXZf3c/N2/K4ZtWrWTzYsVkocqrV8NfasWKCGdNAvDmDbx4Ecmsn8BA6N1bVr+sXVvu0Zsq1efT7vfdKTavGKefnGZZo2XMrD0TUxPTuH/25EgIkWBfJUqUEIqiGMC7d+JB8UZCgDhRsJbIPcZFMBph3egXMX+RX5imc+cKkTOnEAcPCqHTCfHmjRBjxgiRJYsQ3t5hL+vhIUTlykKkSCFEqlRCFC8uxM6doRo8fixEuXJCgBD9+wsRFPT51Hv/96LXjl5CG62JPLPyiCveV+Lt4yclwBkRSU41elKP6qUSvqIYUEiIeNZ7vAhBE09ssouOLRsIRiOyT88uNl/fLHQ63eemK1YIkT+/EHZ2QtjYCNGihRD374e93PnzQjg5CbFsmRCBgUKEhIj/t3f3sVXVZwDHvw+Fljd5EQoiUMCuUF6cRopUnbxsUF7+GBIYdizDkEWDwzH/MBMkGSTIBvjCRsh0hFQystEtoFNeFJkGmGVMKCLWEt4FSse6CXUIaNPeZ3/8LtjAbXuAc8/tvef5JE167jk993l6m6e/+7u/8xzduFH1jjtU394SUV23TrVbN9V27VSLi6/+XCQS0c2HN2vW8iyVhaJztszRC19fCOq3kPSs4BtjvNuxQzUnRxX0zLQJ+uBSN9ovWFug+yr3XT3syuj+8uXYp5kyRXXFiusff3fVCS3pON6Vn2HDVMu+GbnvPr1bR68ZrSxEc1fmasmpEr+zS3mNFXxx+5unvLw83bt3b6LDMCZ8Ll92V7i++CLaujWlk4dT2GcPx1pdoCC7gLkPzWVU31GNNijr0MFdzNX1yj3Djx+HZcvQ117jy5p0WvxqMe1+MZtIC2HbsW2s3LOSTYc3kdk2k/kPz2dW3iwyWmYEk28KEZFSVc2Lue9WCr6I/ABYCAwE7lfVmNVZRMYDvwXSgNWqusTL+a3gG5NgBw/C889DcTGakUHZiFyW3XmCP/espl/3/jw6+FEKhxQyKHPQdT/apQsc2PUlPcu2uqu11q+Hli2pmzGTgX98jqID59heuYmij4o4UX2CzLaZzBk+h6fzn6Z9egM34jVNimfBHwhEgN8Dz8Qq+CKSBhwGxgIVwB7gh6pa3tT5reAb00wcOuQasK1fD9XVfN2+DQf6tmFnh3OUZUJa59v5VuYA+ncZwF2X0uleUU3V1koGnf2QVnU11HbuyPFJI3nn+wP505HTfPTF+9SknwVgZJ+RzMqbxeTcyTai90HcCn69J9hOwwX/AWChqo6Lbs8DUNVfN3VeK/jGNO78edi5E9LSYPTomG3i/VVTA++959Zq7tuHlpcjX3113WEXW8HhLrCjD7wxEEp6Q110JWWLC70ZM+Ahpg8bT0F2gbVF8FljBT+Ie9/0BE7X264AYlyn54jIE8ATAFlZWfGNzJgkpQqLFsHy5ZCf73rZPPaYm3158sk4PnF6ulsrH11wL3V1bm7+4kWoraX60jkq2kU4fVuEsxf/Te3nQqtNrWn/bmsi57MYl9efBfPaM2RIHGM0DWqy4IvI34BYzSrmq+qbHp4j1qc6Db6tUNVVwCpwI3wP5zcmdF55xXW0LC+HHtEB8pEjrtFZjx7wyCMBBZKWBjk5Vzc7Rb/q1/Nnxng7VVkZrFkDVVVw772uQdvt1u3YV01eaauqY1R1SIwvL8Ue3Ii+d73tXkDlzQRrjHE3gHrhBdeJoEe92ZCcHHj5Zbcv2bz0kuvZ06aNa7u8f7/rnVZamujIUksQUzp7gBwR6QecAQqB6QE8rzEpqbrazd0PG3b9voIC18Uymezb5/5RlZa67s3gRvcbNsC0ae6dSwtrAuOLW/o1ishkEakAHgA2i8jW6ON3isgWAFWtBZ4CtgIHgb+oaoxOHMYYL9q1c03KYt1w5LPP6q17TxJFRTB79jfF/oopU1yzzO3bExJWSrqlgq+qb6hqL1XNUNXuV1biqGqlqk6sd9wWVe2vqtmquvhWgzYmzDIyYOpUWHLN1SyqsHix+/A2mZw5AwMGxN6XmwuVNgHsG3ujZEwSWroUNm929wnZsAHWrYOxY92CmXnzEh3djRkyBD744PrH6+qgpAQGDw4+plQVxBy+McZn3bq5m42vXetWtqSluTsCFha6+8wmk8cfh6FD3T+vESPcY5EILFgAffu6lszGH9ZLxxiTcNu2uX9Yd98N2dnu/iddu7qlp927Jzq65JLoC6+MMaZRY8fCyZNumqqqyq3Syc+HRnqzmZtgBd8Y0yy0bu1W5pj4sQ9tjTEmJKzgG2NMSFjBN8aYkLCCb4wxIWEF3xhjQsIKvjHGhESzvvBKRP4DnPThVF2BGK2mUlaY8g1TrhCufMOUK/iXbx9VzYy1o1kXfL+IyN6GrjxLRWHKN0y5QrjyDVOuEEy+NqVjjDEhYQXfGGNCIiwFf1WiAwhYmPINU64QrnzDlCsEkG8o5vCNMcaEZ4RvjDGhZwXfGGNCIqUKvoiMF5FDInJURObG2C8isiK6/4CI3JeIOP3gIdcfRXM8ICK7ROSeRMTpl6byrXfcMBGpE5GpQcbnJy+5isgoEdkvIp+KyI6gY/STh7/ljiKyUUQ+juY7MxFx+kFEikSkSkTKGtgf3xqlqinxBaQBx4C7gHTgY2DQNcdMBN4GBMgH/pnouOOY64NA5+j3E5I1V6/51jvufWALMDXRccfxte0ElANZ0e1uiY47zvk+ByyNfp8JnAPSEx37TeY7ArgPKGtgf1xrVCqN8O8HjqrqcVWtAYqBSdccMwn4gzq7gU4i0iPoQH3QZK6quktVz0c3dwO9Ao7RT15eW4CfARuAqiCD85mXXKcDr6vqKQBVTfV8FbhNRARojyv4tcGG6Q9V3YmLvyFxrVGpVPB7AqfrbVdEH7vRY5LBjebxE9yoIVk1ma+I9AQmA68GGFc8eHlt+wOdRWS7iJSKyIzAovOfl3xXAgOBSuAT4OeqGgkmvMDFtUal0i0OY9398to1p16OSQae8xCR0biC/524RhRfXvL9DfCsqtZJct8I1UuuLYGhwPeANsA/RGS3qh6Od3Bx4CXfccB+4LtANrBNRP6uqv+Ld3AJENcalUoFvwLoXW+7F25EcKPHJANPeYjIt4HVwARV/Tyg2OLBS755QHG02HcFJopIrar+NZgQfeP17/i/qnoRuCgiO4F7gGQs+F7ynQksUTfJfVRETgC5wIfBhBiouNaoVJrS2QPkiEg/EUkHCoG3rjnmLWBG9JPwfOALVf1X0IH6oMlcRSQLeB34cZKO/OprMl9V7aeqfVW1L7Ae+GkSFnvw9nf8JvCwiLQUkbbAcOBgwHH6xUu+p3DvZhCR7sAA4HigUQYnrjUqZUb4qlorIk8BW3Gf/Bep6qciMiu6/1Xc6o2JwFHgEm7kkHQ85vpLoAvwu+iot1aTtPOgx3xTgpdcVfWgiLwDHAAiwGpVjbnMr7nz+NouAtaIyCe4KY9nVTUp2yaLyDpgFNBVRCqABUArCKZGWWsFY4wJiVSa0jHGGNMIK/jGGBMSVvCNMSYkrOAbY0xIWME3xpiQsIJvjDEhYQXfGGNC4v/oM7oB7mv66gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate data\n",
    "feature = PloynomialFeature(7)\n",
    "X_train = feature.transform(x_train)\n",
    "X_test = feature.transform(x_test)\n",
    "\n",
    "# Model\n",
    "model = LinearRegression()\n",
    "\n",
    "# training\n",
    "model.fit(X_train, y_train)\n",
    "\n",
    "# testing\n",
    "y = model.predict(X_test)\n",
    "\n",
    "# plot\n",
    "plt.scatter(x_train, y_train, fc='none', ec='b', s=50, label='train data')\n",
    "plt.plot(x_exact, y_exact, c='g', label='$\\sin 2\\pi x$')\n",
    "plt.plot(x_test, y, c='r', label='fitting')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1,2])\n",
    "print(x.ndim)\n",
    "print(x[:,None].ndim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Regression\n",
    "$Ols \\longrightarrow \\omega_{MLE} \\stackrel{x}{\\longrightarrow} y$\n",
    "\n",
    "$Ridge \\longrightarrow \\omega_{MAP} \\stackrel{x}{\\longrightarrow} y$\n",
    "\n",
    "$Bayes \\longrightarrow p(\\omega|D) \\stackrel{x}{\\longrightarrow} p(y|x,D)(代表着y的范围)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BayesianRegression(Regression):\n",
    "    def __init__(self, a=1., b=1.):\n",
    "        self.a = a\n",
    "        self.b = b\n",
    "        self.w_mean = None # mu_0\n",
    "        self.w_prec = None # Lambda_0\n",
    "        \n",
    "    def _is_prior_defined(self):\n",
    "        return self.w_mean is not None and self.w_prec is not None\n",
    "        # 都存在，认为先验分布已经定义了\n",
    "        \n",
    "    def _get_prior(self):\n",
    "        if self._is_prior_defined():\n",
    "            return self.w_mean, self.w_prec\n",
    "        else:\n",
    "            return np.zeros(ndim), 1 / self.b * np.eyes(ndim)\n",
    "        \n",
    "    def fit(self):\n",
    "        \"\"\"\n",
    "        params\n",
    "        ======\n",
    "        X: ndarray with shape(N, K)\n",
    "        \"\"\"\n",
    "        mean_prior, prec_prior = self._get_prior(self, X.shape[1])\n",
    "        \n",
    "        w_prec = prec_prior + self.a * X.T @ X\n",
    "        \n",
    "        w_mean = np.linalg.solve(w_prec, prec_prior @ mean_prior + self.a * X.T @ y)\n",
    "        # 希望数据一个一个进，如果是一次性可以按照课件写\n",
    "        self.w_mean = w_mean\n",
    "        self.w_prec = w_prec\n",
    "        self.W_cov = np.linalg.inv(self.w_prec)\n",
    "        \n",
    "    def predict(self,X):\n",
    "        return X @ self.w_mean\n",
    "    \n",
    "    \n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.69106813,  1.29031405,  1.81605884],\n",
       "       [-0.95460617,  0.67594761,  0.63450562],\n",
       "       [-1.32161986,  1.10402855, -2.53273872]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.random.randn(3,3)\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.69106813, -0.95460617, -1.32161986],\n",
       "       [ 1.29031405,  0.67594761,  1.10402855],\n",
       "       [ 1.81605884,  0.63450562, -2.53273872]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.transpose()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
